home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Gold Collection / Software Vault - The Gold Collection (American Databankers) (1993).ISO / cdr53 / acmalg01.zip / ACM691.FOR < prev    next >
Text File  |  1993-01-01  |  185KB  |  4,692 lines
  1. C      ALGORITHM 691, COLLECTED ALGORITHMS FROM ACM.
  2. C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
  3. C      VOL. 17, NO. 2, PP. 218-232.  JUNE, 1991.
  4. This disk contains 8 files:
  5.  
  6. README : this information file
  7.  
  8. stest.for : which contains simple test routines (single precision)
  9.  
  10.           TEST    (main)
  11.           F
  12.  
  13. squad.for : which contains quadrature routines (single precision)
  14.  
  15.           QXGS
  16.           QXGSE
  17.           QXG
  18.           QXGE
  19.           QXLQM
  20.           QXRUL
  21.           QXRRD
  22.           QXCPY
  23.  
  24. saux.for : which contains auxiliary routines (single precision)
  25.  
  26.           QPSRT
  27.           QELG
  28.           R1MACH (machine dependent, TO BE MODIFIED)
  29.  
  30. dtest.for : which contains simple test routines (double precision)
  31.  
  32.           DTEST    (main)
  33.           F
  34.  
  35. dquad.for : which contains quadrature routines (double precision)
  36.  
  37.           DQXGS
  38.           DQXGSE
  39.           DQXG
  40.           DQXGE
  41.           DQXLQM
  42.           DQXRUL
  43.           DQXRRD
  44.           DQXCPY
  45.  
  46. daux.for : which contains auxiliary routines (double precision)
  47.  
  48.           DQPSRT
  49.           DQELG
  50.           D1MACH  (machine dependent, TO BE MODIFIED)
  51.  
  52. CC+------------------------------------------------------------------+CC
  53. CC+                                                                  +CC
  54. CC+                LAWRENCE LIVERMORE NATIONAL LABORATORY            +CC
  55. CC+                                                                  +CC
  56. CC+        MATHEMATICS AND STATISTICS DIVISION SOFTWARE LIBRARY      +CC
  57. CC+                                                                  +CC
  58. CC+------------------------------------------------------------------+CC
  59. CC+                                                                  +CC
  60. CC+  THE SLATEC MATHEMATICAL PROGRAM LIBRARY IS ISSUED BY            +CC
  61. CC+  THE FOLLOWING:                                                  +CC
  62. CC+                                                                  +CC
  63. CC+       AIR FORCE WEAPONS LABORATORY, ALBUQUERQUE                  +CC
  64. CC+       LAWRENCE LIVERMORE NATIONAL LABORATORY, LIVERMORE          +CC
  65. CC+       LOS ALAMOS NATIONAL SCIENTIFIC LABORATORY, LOS ALAMOS      +CC
  66. CC+       NATIONAL BUREAU OF STANDARDS, WASHINGTON                   +CC
  67. CC+       OAK RIDGE NATIONAL LABORATORY, OAK RIDGE                   +CC
  68. CC+       SANDIA NATIONAL LABORATORIES, ALBUQUERQUE                  +CC
  69. CC+       SANDIA NATIONAL LABORATORIES, LIVERMORE                    +CC
  70. CC+                                                                  +CC
  71. CC+  PLEASE REFER QUESTIONS REGARDING THIS SOFTWARE TO THE           +CC
  72. CC+  MSD CONSULTING OFFICE, EXTENSION 3-2976.                        +CC
  73. CC+                                                                  +CC
  74. CC+                                                                  +CC
  75. CC+              * * * * * N O T I C E * * * * *                     +CC
  76. CC+                                                                  +CC
  77. CC+  THIS MATERIAL WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED      +CC
  78. CC+  BY THE UNITED STATES GOVERNMENT.  NEITHER THE UNITED            +CC
  79. CC+  STATES GOVERNMENT, NOR THE DEPARTMENT OF DEFENSE, NOR           +CC
  80. CC+  ANY OF THEIR EMPLOYEES, MAKES ANY WARRANTY, EXPRESSED           +CC
  81. CC+  OR IMPLIED, OR ASSUMES ANY LEGAL LIABILITY OR                   +CC
  82. CC+  RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR               +CC
  83. CC+  USEFULNESS OF ANY INFORMATION, APPARATUS, PRODUCT,              +CC
  84. CC+  OR PROCESS DISCLOSED, OR REPRESENTS ITS USE WOULD NOT           +CC
  85. CC+  INFRINGE UPON PRIVATELY OWNED RIGHTS.                           +CC
  86. CC+                                                                  +CC
  87. CC+------------------------------------------------------------------+CC
  88. *DECK DQPSRT
  89.       SUBROUTINE DQPSRT(LIMIT,LAST,MAXERR,ERMAX,ELIST,IORD,NRMAX)
  90. C***BEGIN PROLOGUE  DQPSRT
  91. C***REFER TO  DQAGE,DQAGIE,DQAGPE,DQAWSE
  92. C***ROUTINES CALLED  (NONE)
  93. C***REVISION DATE  810101   (YYMMDD)
  94. C***KEYWORDS  SEQUENTIAL SORTING
  95. C***AUTHOR  PIESSENS, ROBERT, APPLIED MATH. AND PROGR. DIV. -
  96. C             K. U. LEUVEN
  97. C           DE DONCKER, ELISE, APPLIED MATH. AND PROGR. DIV. -
  98. C             K. U. LEUVEN
  99. C***PURPOSE  THIS ROUTINE MAINTAINS THE DESCENDING ORDERING IN THE
  100. C            LIST OF THE LOCAL ERROR ESTIMATED RESULTING FROM THE
  101. C            INTERVAL SUBDIVISION PROCESS. AT EACH CALL TWO ERROR
  102. C            ESTIMATES ARE INSERTED USING THE SEQUENTIAL SEARCH
  103. C            METHOD, TOP-DOWN FOR THE LARGEST ERROR ESTIMATE AND
  104. C            BOTTOM-UP FOR THE SMALLEST ERROR ESTIMATE.
  105. C***DESCRIPTION
  106. C
  107. C           ORDERING ROUTINE
  108. C           STANDARD FORTRAN SUBROUTINE
  109. C           DOUBLE PRECISION VERSION
  110. C
  111. C           PARAMETERS (MEANING AT OUTPUT)
  112. C              LIMIT  - INTEGER
  113. C                       MAXIMUM NUMBER OF ERROR ESTIMATES THE LIST
  114. C                       CAN CONTAIN
  115. C
  116. C              LAST   - INTEGER
  117. C                       NUMBER OF ERROR ESTIMATES CURRENTLY IN THE LIST
  118. C
  119. C              MAXERR - INTEGER
  120. C                       MAXERR POINTS TO THE NRMAX-TH LARGEST ERROR
  121. C                       ESTIMATE CURRENTLY IN THE LIST
  122. C
  123. C              ERMAX  - DOUBLE PRECISION
  124. C                       NRMAX-TH LARGEST ERROR ESTIMATE
  125. C                       ERMAX = ELIST(MAXERR)
  126. C
  127. C              ELIST  - DOUBLE PRECISION
  128. C                       VECTOR OF DIMENSION LAST CONTAINING
  129. C                       THE ERROR ESTIMATES
  130. C
  131. C              IORD   - INTEGER
  132. C                       VECTOR OF DIMENSION LAST, THE FIRST K ELEMENTS
  133. C                       OF WHICH CONTAIN POINTERS TO THE ERROR
  134. C                       ESTIMATES, SUCH THAT
  135. C                       ELIST(IORD(1)),...,  ELIST(IORD(K))
  136. C                       FORM A DECREASING SEQUENCE, WITH
  137. C                       K = LAST IF LAST.LE.(LIMIT/2+2), AND
  138. C                       K = LIMIT+1-LAST OTHERWISE
  139. C
  140. C              NRMAX  - INTEGER
  141. C                       MAXERR = IORD(NRMAX)
  142. C***END PROLOGUE  DQPSRT
  143. C
  144.       DOUBLE PRECISION ELIST,ERMAX,ERRMAX,ERRMIN
  145.       INTEGER I,IBEG,IDO,IORD,ISUCC,J,JBND,JUPBN,K,LAST,LIMIT,MAXERR,
  146.      1  NRMAX
  147.       DIMENSION ELIST(LAST),IORD(LAST)
  148. C
  149. C           CHECK WHETHER THE LIST CONTAINS MORE THAN
  150. C           TWO ERROR ESTIMATES.
  151. C
  152. C***FIRST EXECUTABLE STATEMENT  DQPSRT
  153.       IF(LAST.GT.2) GO TO 10
  154.       IORD(1) = 1
  155.       IORD(2) = 2
  156.       GO TO 90
  157. C
  158. C           THIS PART OF THE ROUTINE IS ONLY EXECUTED IF, DUE TO A
  159. C           DIFFICULT INTEGRAND, SUBDIVISION INCREASED THE ERROR
  160. C           ESTIMATE. IN THE NORMAL CASE THE INSERT PROCEDURE SHOULD
  161. C           START AFTER THE NRMAX-TH LARGEST ERROR ESTIMATE.
  162. C
  163.    10 ERRMAX = ELIST(MAXERR)
  164.       IF(NRMAX.EQ.1) GO TO 30
  165.       IDO = NRMAX-1
  166.       DO 20 I = 1,IDO
  167.         ISUCC = IORD(NRMAX-1)
  168. C ***JUMP OUT OF DO-LOOP
  169.         IF(ERRMAX.LE.ELIST(ISUCC)) GO TO 30
  170.         IORD(NRMAX) = ISUCC
  171.         NRMAX = NRMAX-1
  172.    20    CONTINUE
  173. C
  174. C           COMPUTE THE NUMBER OF ELEMENTS IN THE LIST TO BE MAINTAINED
  175. C           IN DESCENDING ORDER. THIS NUMBER DEPENDS ON THE NUMBER OF
  176. C           SUBDIVISIONS STILL ALLOWED.
  177. C
  178.    30 JUPBN = LAST
  179.       IF(LAST.GT.(LIMIT/2+2)) JUPBN = LIMIT+3-LAST
  180.       ERRMIN = ELIST(LAST)
  181. C
  182. C           INSERT ERRMAX BY TRAVERSING THE LIST TOP-DOWN,
  183. C           STARTING COMPARISON FROM THE ELEMENT ELIST(IORD(NRMAX+1)).
  184. C
  185.       JBND = JUPBN-1
  186.       IBEG = NRMAX+1
  187.       IF(IBEG.GT.JBND) GO TO 50
  188.       DO 40 I=IBEG,JBND
  189.         ISUCC = IORD(I)
  190. C ***JUMP OUT OF DO-LOOP
  191.         IF(ERRMAX.GE.ELIST(ISUCC)) GO TO 60
  192.         IORD(I-1) = ISUCC
  193.    40 CONTINUE
  194.    50 IORD(JBND) = MAXERR
  195.       IORD(JUPBN) = LAST
  196.       GO TO 90
  197. C
  198. C           INSERT ERRMIN BY TRAVERSING THE LIST BOTTOM-UP.
  199. C
  200.    60 IORD(I-1) = MAXERR
  201.       K = JBND
  202.       DO 70 J=I,JBND
  203.         ISUCC = IORD(K)
  204. C ***JUMP OUT OF DO-LOOP
  205.         IF(ERRMIN.LT.ELIST(ISUCC)) GO TO 80
  206.         IORD(K+1) = ISUCC
  207.         K = K-1
  208.    70 CONTINUE
  209.       IORD(I) = LAST
  210.       GO TO 90
  211.    80 IORD(K+1) = LAST
  212. C
  213. C           SET MAXERR AND ERMAX.
  214. C
  215.    90 MAXERR = IORD(NRMAX)
  216.       ERMAX = ELIST(MAXERR)
  217.       RETURN
  218.       END
  219. *DECK DQELG
  220.       SUBROUTINE DQELG(N,EPSTAB,RESULT,ABSERR,RES3LA,NRES)
  221. C***BEGIN PROLOGUE  DQELG
  222. C***REFER TO  DQAGIE,DQAGOE,DQAGPE,DQAGSE
  223. C***ROUTINES CALLED  D1MACH
  224. C***REVISION DATE  830518   (YYMMDD)
  225. C***KEYWORDS  CONVERGENCE ACCELERATION,EPSILON ALGORITHM,EXTRAPOLATION
  226. C***AUTHOR  PIESSENS, ROBERT, APPLIED MATH. AND PROGR. DIV. -
  227. C             K. U. LEUVEN
  228. C           DE DONCKER, ELISE, APPLIED MATH. AND PROGR. DIV. -
  229. C             K. U. LEUVEN
  230. C***PURPOSE  THE ROUTINE DETERMINES THE LIMIT OF A GIVEN SEQUENCE OF
  231. C            APPROXIMATIONS, BY MEANS OF THE EPSILON ALGORITHM OF
  232. C            P.WYNN. AN ESTIMATE OF THE ABSOLUTE ERROR IS ALSO GIVEN.
  233. C            THE CONDENSED EPSILON TABLE IS COMPUTED. ONLY THOSE
  234. C            ELEMENTS NEEDED FOR THE COMPUTATION OF THE NEXT DIAGONAL
  235. C            ARE PRESERVED.
  236. C***DESCRIPTION
  237. C
  238. C           EPSILON ALGORITHM
  239. C           STANDARD FORTRAN SUBROUTINE
  240. C           DOUBLE PRECISION VERSION
  241. C
  242. C           PARAMETERS
  243. C              N      - INTEGER
  244. C                       EPSTAB(N) CONTAINS THE NEW ELEMENT IN THE
  245. C                       FIRST COLUMN OF THE EPSILON TABLE.
  246. C
  247. C              EPSTAB - DOUBLE PRECISION
  248. C                       VECTOR OF DIMENSION 52 CONTAINING THE ELEMENTS
  249. C                       OF THE TWO LOWER DIAGONALS OF THE TRIANGULAR
  250. C                       EPSILON TABLE. THE ELEMENTS ARE NUMBERED
  251. C                       STARTING AT THE RIGHT-HAND CORNER OF THE
  252. C                       TRIANGLE.
  253. C
  254. C              RESULT - DOUBLE PRECISION
  255. C                       RESULTING APPROXIMATION TO THE INTEGRAL
  256. C
  257. C              ABSERR - DOUBLE PRECISION
  258. C                       ESTIMATE OF THE ABSOLUTE ERROR COMPUTED FROM
  259. C                       RESULT AND THE 3 PREVIOUS RESULTS
  260. C
  261. C              RES3LA - DOUBLE PRECISION
  262. C                       VECTOR OF DIMENSION 3 CONTAINING THE LAST 3
  263. C                       RESULTS
  264. C
  265. C              NRES   - INTEGER
  266. C                       NUMBER OF CALLS TO THE ROUTINE
  267. C                       (SHOULD BE ZERO AT FIRST CALL)
  268. C***END PROLOGUE  DQELG
  269. C
  270.       DOUBLE PRECISION ABSERR,DABS,DELTA1,DELTA2,DELTA3,DMAX1,D1MACH,
  271.      1  EPMACH,EPSINF,EPSTAB,ERROR,ERR1,ERR2,ERR3,E0,E1,E1ABS,E2,E3,
  272.      2  OFLOW,RES,RESULT,RES3LA,SS,TOL1,TOL2,TOL3
  273.       INTEGER I,IB,IB2,IE,INDX,K1,K2,K3,LIMEXP,N,NEWELM,NRES,NUM
  274.       DIMENSION EPSTAB(52),RES3LA(3)
  275. C
  276. C           LIST OF MAJOR VARIABLES
  277. C           -----------------------
  278. C
  279. C           E0     - THE 4 ELEMENTS ON WHICH THE COMPUTATION OF A NEW
  280. C           E1       ELEMENT IN THE EPSILON TABLE IS BASED
  281. C           E2
  282. C           E3                 E0
  283. C                        E3    E1    NEW
  284. C                              E2
  285. C           NEWELM - NUMBER OF ELEMENTS TO BE COMPUTED IN THE NEW
  286. C                    DIAGONAL
  287. C           ERROR  - ERROR = ABS(E1-E0)+ABS(E2-E1)+ABS(NEW-E2)
  288. C           RESULT - THE ELEMENT IN THE NEW DIAGONAL WITH LEAST VALUE
  289. C                    OF ERROR
  290. C
  291. C           MACHINE DEPENDENT CONSTANTS
  292. C           ---------------------------
  293. C
  294. C           EPMACH IS THE LARGEST RELATIVE SPACING.
  295. C           OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
  296. C           LIMEXP IS THE MAXIMUM NUMBER OF ELEMENTS THE EPSILON
  297. C           TABLE CAN CONTAIN. IF THIS NUMBER IS REACHED, THE UPPER
  298. C           DIAGONAL OF THE EPSILON TABLE IS DELETED.
  299. C
  300. C***FIRST EXECUTABLE STATEMENT  DQELG
  301.       EPMACH = D1MACH(4)
  302.       OFLOW = D1MACH(2)
  303.       NRES = NRES+1
  304.       ABSERR = OFLOW
  305.       RESULT = EPSTAB(N)
  306.       IF(N.LT.3) GO TO 100
  307.       LIMEXP = 50
  308.       EPSTAB(N+2) = EPSTAB(N)
  309.       NEWELM = (N-1)/2
  310.       EPSTAB(N) = OFLOW
  311.       NUM = N
  312.       K1 = N
  313.       DO 40 I = 1,NEWELM
  314.         K2 = K1-1
  315.         K3 = K1-2
  316.         RES = EPSTAB(K1+2)
  317.         E0 = EPSTAB(K3)
  318.         E1 = EPSTAB(K2)
  319.         E2 = RES
  320.         E1ABS = DABS(E1)
  321.         DELTA2 = E2-E1
  322.         ERR2 = DABS(DELTA2)
  323.         TOL2 = DMAX1(DABS(E2),E1ABS)*EPMACH
  324.         DELTA3 = E1-E0
  325.         ERR3 = DABS(DELTA3)
  326.         TOL3 = DMAX1(E1ABS,DABS(E0))*EPMACH
  327.         IF(ERR2.GT.TOL2.OR.ERR3.GT.TOL3) GO TO 10
  328. C
  329. C           IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE
  330. C           ACCURACY, CONVERGENCE IS ASSUMED.
  331. C           RESULT = E2
  332. C           ABSERR = ABS(E1-E0)+ABS(E2-E1)
  333. C
  334.         RESULT = RES
  335.         ABSERR = ERR2+ERR3
  336. C ***JUMP OUT OF DO-LOOP
  337.         GO TO 100
  338.    10   E3 = EPSTAB(K1)
  339.         EPSTAB(K1) = E1
  340.         DELTA1 = E1-E3
  341.         ERR1 = DABS(DELTA1)
  342.         TOL1 = DMAX1(E1ABS,DABS(E3))*EPMACH
  343. C
  344. C           IF TWO ELEMENTS ARE VERY CLOSE TO EACH OTHER, OMIT
  345. C           A PART OF THE TABLE BY ADJUSTING THE VALUE OF N
  346. C
  347.         IF(ERR1.LE.TOL1.OR.ERR2.LE.TOL2.OR.ERR3.LE.TOL3) GO TO 20
  348.         SS = 0.1D+01/DELTA1+0.1D+01/DELTA2-0.1D+01/DELTA3
  349.         EPSINF = DABS(SS*E1)
  350. C
  351. C           TEST TO DETECT IRREGULAR BEHAVIOUR IN THE TABLE, AND
  352. C           EVENTUALLY OMIT A PART OF THE TABLE ADJUSTING THE VALUE
  353. C           OF N.
  354. C
  355.         IF(EPSINF.GT.0.1D-03) GO TO 30
  356.    20   N = I+I-1
  357. C ***JUMP OUT OF DO-LOOP
  358.         GO TO 50
  359. C
  360. C           COMPUTE A NEW ELEMENT AND EVENTUALLY ADJUST
  361. C           THE VALUE OF RESULT.
  362. C
  363.    30   RES = E1+0.1D+01/SS
  364.         EPSTAB(K1) = RES
  365.         K1 = K1-2
  366.         ERROR = ERR2+DABS(RES-E2)+ERR3
  367.         IF(ERROR.GT.ABSERR) GO TO 40
  368.         ABSERR = ERROR
  369.         RESULT = RES
  370.    40 CONTINUE
  371. C
  372. C           SHIFT THE TABLE.
  373. C
  374.    50 IF(N.EQ.LIMEXP) N = 2*(LIMEXP/2)-1
  375.       IB = 1
  376.       IF((NUM/2)*2.EQ.NUM) IB = 2
  377.       IE = NEWELM+1
  378.       DO 60 I=1,IE
  379.         IB2 = IB+2
  380.         EPSTAB(IB) = EPSTAB(IB2)
  381.         IB = IB2
  382.    60 CONTINUE
  383.       IF(NUM.EQ.N) GO TO 80
  384.       INDX = NUM-N+1
  385.       DO 70 I = 1,N
  386.         EPSTAB(I)= EPSTAB(INDX)
  387.         INDX = INDX+1
  388.    70 CONTINUE
  389.    80 IF(NRES.GE.4) GO TO 90
  390.       RES3LA(NRES) = RESULT
  391.       ABSERR = OFLOW
  392.       GO TO 100
  393. C
  394. C           COMPUTE ERROR ESTIMATE
  395. C
  396.    90 ABSERR = DABS(RESULT-RES3LA(3))+DABS(RESULT-RES3LA(2))
  397.      1  +DABS(RESULT-RES3LA(1))
  398.       RES3LA(1) = RES3LA(2)
  399.       RES3LA(2) = RES3LA(3)
  400.       RES3LA(3) = RESULT
  401.   100 ABSERR = DMAX1(ABSERR,0.5D+01*EPMACH*DABS(RESULT))
  402.       RETURN
  403.       END
  404. *DECK D1MACH
  405.       DOUBLE PRECISION FUNCTION D1MACH(I)
  406. C***BEGIN PROLOGUE  D1MACH
  407. C***DATE WRITTEN   750101   (YYMMDD)
  408. C***REVISION DATE  890213   (YYMMDD)
  409. C***CATEGORY NO.  R1
  410. C***KEYWORDS  LIBRARY=SLATEC,TYPE=DOUBLE PRECISION(R1MACH-S D1MACH-D),
  411. C             MACHINE CONSTANTS
  412. C***AUTHOR  FOX, P. A., (BELL LABS)
  413. C           HALL, A. D., (BELL LABS)
  414. C           SCHRYER, N. L., (BELL LABS)
  415. C***PURPOSE  RETURNS DOUBLE PRECISION MACHINE DEPENDENT CONSTANTS
  416. C***DESCRIPTION
  417. C
  418. C   D1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS
  419. C   FOR THE LOCAL MACHINE ENVIRONMENT.  IT IS A FUNCTION
  420. C   SUBPROGRAM WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED
  421. C   AS FOLLOWS, FOR EXAMPLE
  422. C
  423. C        D = D1MACH(I)
  424. C
  425. C   WHERE I=1,...,5.  THE (OUTPUT) VALUE OF D ABOVE IS
  426. C   DETERMINED BY THE (INPUT) VALUE OF I.  THE RESULTS FOR
  427. C   VARIOUS VALUES OF I ARE DISCUSSED BELOW.
  428. C
  429. C   D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
  430. C   D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
  431. C   D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
  432. C   D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
  433. C   D1MACH( 5) = LOG10(B)
  434. C
  435. C   ASSUME DOUBLE PRECISION NUMBERS ARE REPRESENTED IN THE T-DIGIT,
  436. C   BASE-B FORM
  437. C
  438. C              SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
  439. C
  440. C   WHERE 0 .LE. X(I) .LT. B FOR I=1,...,T, 0 .LT. X(1), AND
  441. C   EMIN .LE. E .LE. EMAX.
  442. C
  443. C   THE VALUES OF B, T, EMIN AND EMAX ARE PROVIDED IN I1MACH AS
  444. C   FOLLOWS:
  445. C   I1MACH(10) = B, THE BASE.
  446. C   I1MACH(14) = T, THE NUMBER OF BASE-B DIGITS.
  447. C   I1MACH(15) = EMIN, THE SMALLEST EXPONENT E.
  448. C   I1MACH(16) = EMAX, THE LARGEST EXPONENT E.
  449. C
  450. C   TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT,
  451. C   THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY
  452. C   REMOVING THE C FROM COLUMN 1.  ALSO, THE VALUES OF
  453. C   D1MACH(1) - D1MACH(4) SHOULD BE CHECKED FOR CONSISTENCY
  454. C   WITH THE LOCAL OPERATING SYSTEM.
  455. C
  456. C***REFERENCES  FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A
  457. C                 PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL
  458. C                 SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188.
  459. C***ROUTINES CALLED   XERROR
  460. C
  461. C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
  462. C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
  463. C        NOTE ADDED BY F. ROMANI  7/11/89
  464. C
  465. C***END PROLOGUE  D1MACH
  466. C
  467.       INTEGER SMALL(4)
  468.       INTEGER LARGE(4)
  469.       INTEGER RIGHT(4)
  470.       INTEGER DIVER(4)
  471.       INTEGER LOG10(4)
  472. C
  473.       DOUBLE PRECISION DMACH(5)
  474.       SAVE DMACH
  475. C
  476.       EQUIVALENCE (DMACH(1),SMALL(1))
  477.       EQUIVALENCE (DMACH(2),LARGE(1))
  478.       EQUIVALENCE (DMACH(3),RIGHT(1))
  479.       EQUIVALENCE (DMACH(4),DIVER(1))
  480.       EQUIVALENCE (DMACH(5),LOG10(1))
  481. C
  482. C     MACHINE CONSTANTS FOR THE AMIGA
  483. C     ABSOFT FORTRAN COMPILER USING THE 68020/68881 COMPILER OPTION
  484. C
  485. C     DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' /
  486. C     DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' /
  487. C     DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' /
  488. C     DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' /
  489. C     DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' /
  490. C
  491. C     MACHINE CONSTANTS FOR THE AMIGA
  492. C     ABSOFT FORTRAN COMPILER USING SOFTWARE FLOATING POINT
  493. C
  494. C     DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' /
  495. C     DATA LARGE(1), LARGE(2) / Z'7FDFFFFF', Z'FFFFFFFF' /
  496. C     DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' /
  497. C     DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' /
  498. C     DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' /
  499. C
  500. C     MACHINE CONSTANTS FOR THE APOLLO
  501. C
  502. C     DATA SMALL(1), SMALL(2) / 16#00100000, 16#00000000 /
  503. C     DATA LARGE(1), LARGE(2) / 16#7FFFFFFF, 16#FFFFFFFF /
  504. C     DATA RIGHT(1), RIGHT(2) / 16#3CA00000, 16#00000000 /
  505. C     DATA DIVER(1), DIVER(2) / 16#3CB00000, 16#00000000 /
  506. C     DATA LOG10(1), LOG10(2) / 16#3FD34413, 16#509F79FF /
  507. C
  508. C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM
  509. C
  510. C     DATA SMALL(1) / ZC00800000 /
  511. C     DATA SMALL(2) / Z000000000 /
  512. C     DATA LARGE(1) / ZDFFFFFFFF /
  513. C     DATA LARGE(2) / ZFFFFFFFFF /
  514. C     DATA RIGHT(1) / ZCC5800000 /
  515. C     DATA RIGHT(2) / Z000000000 /
  516. C     DATA DIVER(1) / ZCC6800000 /
  517. C     DATA DIVER(2) / Z000000000 /
  518. C     DATA LOG10(1) / ZD00E730E7 /
  519. C     DATA LOG10(2) / ZC77800DC0 /
  520. C
  521. C     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM
  522. C
  523. C     DATA SMALL(1) / O1771000000000000 /
  524. C     DATA SMALL(2) / O0000000000000000 /
  525. C     DATA LARGE(1) / O0777777777777777 /
  526. C     DATA LARGE(2) / O0007777777777777 /
  527. C     DATA RIGHT(1) / O1461000000000000 /
  528. C     DATA RIGHT(2) / O0000000000000000 /
  529. C     DATA DIVER(1) / O1451000000000000 /
  530. C     DATA DIVER(2) / O0000000000000000 /
  531. C     DATA LOG10(1) / O1157163034761674 /
  532. C     DATA LOG10(2) / O0006677466732724 /
  533. C
  534. C     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS
  535. C
  536. C     DATA SMALL(1) / O1771000000000000 /
  537. C     DATA SMALL(2) / O7770000000000000 /
  538. C     DATA LARGE(1) / O0777777777777777 /
  539. C     DATA LARGE(2) / O7777777777777777 /
  540. C     DATA RIGHT(1) / O1461000000000000 /
  541. C     DATA RIGHT(2) / O0000000000000000 /
  542. C     DATA DIVER(1) / O1451000000000000 /
  543. C     DATA DIVER(2) / O0000000000000000 /
  544. C     DATA LOG10(1) / O1157163034761674 /
  545. C     DATA LOG10(2) / O0006677466732724 /
  546. C
  547. C     MACHINE CONSTANTS FOR THE CDC 170/180 SERIES USING NOS/VE
  548. C
  549. C     DATA SMALL(1) / Z"3001800000000000" /
  550. C     DATA SMALL(2) / Z"3001000000000000" /
  551. C     DATA LARGE(1) / Z"4FFEFFFFFFFFFFFE" /
  552. C     DATA LARGE(2) / Z"4FFE000000000000" /
  553. C     DATA RIGHT(1) / Z"3FD2800000000000" /
  554. C     DATA RIGHT(2) / Z"3FD2000000000000" /
  555. C     DATA DIVER(1) / Z"3FD3800000000000" /
  556. C     DATA DIVER(2) / Z"3FD3000000000000" /
  557. C     DATA LOG10(1) / Z"3FFF9A209A84FBCF" /
  558. C     DATA LOG10(2) / Z"3FFFF7988F8959AC" /
  559. C
  560. C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES
  561. C
  562. C     DATA SMALL(1) / 00564000000000000000B /
  563. C     DATA SMALL(2) / 00000000000000000000B /
  564. C     DATA LARGE(1) / 37757777777777777777B /
  565. C     DATA LARGE(2) / 37157777777777777777B /
  566. C     DATA RIGHT(1) / 15624000000000000000B /
  567. C     DATA RIGHT(2) / 00000000000000000000B /
  568. C     DATA DIVER(1) / 15634000000000000000B /
  569. C     DATA DIVER(2) / 00000000000000000000B /
  570. C     DATA LOG10(1) / 17164642023241175717B /
  571. C     DATA LOG10(2) / 16367571421742254654B /
  572. C
  573. C     MACHINE CONSTANTS FOR THE CELERITY C1260
  574. C
  575. C     DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' /
  576. C     DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' /
  577. C     DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' /
  578. C     DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' /
  579. C     DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' /
  580. C
  581. C     MACHINE CONSTANTS FOR THE CONVEX C-1
  582. C
  583. C     DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X /
  584. C     DATA LARGE(1), LARGE(2) / '7FFFFFFF'X,'FFFFFFFF'X /
  585. C     DATA RIGHT(1), RIGHT(2) / '3CC00000'X,'00000000'X /
  586. C     DATA DIVER(1), DIVER(2) / '3CD00000'X,'00000000'X /
  587. C     DATA LOG10(1), LOG10(2) / '3FF34413'X,'509F79FF'X /
  588. C
  589. C     MACHINE CONSTANTS FOR THE CRAY-1
  590. C
  591. C     DATA SMALL(1) / 201354000000000000000B /
  592. C     DATA SMALL(2) / 000000000000000000000B /
  593. C     DATA LARGE(1) / 577767777777777777777B /
  594. C     DATA LARGE(2) / 000007777777777777774B /
  595. C     DATA RIGHT(1) / 376434000000000000000B /
  596. C     DATA RIGHT(2) / 000000000000000000000B /
  597. C     DATA DIVER(1) / 376444000000000000000B /
  598. C     DATA DIVER(2) / 000000000000000000000B /
  599. C     DATA LOG10(1) / 377774642023241175717B /
  600. C     DATA LOG10(2) / 000007571421742254654B /
  601. C
  602. C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
  603. C
  604. C     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD -
  605. C     STATIC DMACH(5)
  606. C
  607. C     DATA SMALL /    20K, 3*0 /
  608. C     DATA LARGE / 77777K, 3*177777K /
  609. C     DATA RIGHT / 31420K, 3*0 /
  610. C     DATA DIVER / 32020K, 3*0 /
  611. C     DATA LOG10 / 40423K, 42023K, 50237K, 74776K /
  612. C
  613. C     MACHINE CONSTANTS FOR THE ELXSI 6400
  614. C     (ASSUMING REAL*8 IS THE DEFAULT DOUBLE PRECISION)
  615. C
  616. C     DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X /
  617. C     DATA LARGE(1), LARGE(2) / '7FEFFFFF'X,'FFFFFFFF'X /
  618. C     DATA RIGHT(1), RIGHT(2) / '3CB00000'X,'00000000'X /
  619. C     DATA DIVER(1), DIVER(2) / '3CC00000'X,'00000000'X /
  620. C     DATA LOG10(1), LOG10(2) / '3FD34413'X,'509F79FF'X /
  621. C
  622. C     MACHINE CONSTANTS FOR THE HARRIS 220
  623. C
  624. C     DATA SMALL(1), SMALL(2) / '20000000, '00000201 /
  625. C     DATA LARGE(1), LARGE(2) / '37777777, '37777577 /
  626. C     DATA RIGHT(1), RIGHT(2) / '20000000, '00000333 /
  627. C     DATA DIVER(1), DIVER(2) / '20000000, '00000334 /
  628. C     DATA LOG10(1), LOG10(2) / '23210115, '10237777 /
  629. C
  630. C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES
  631. C
  632. C     DATA SMALL(1), SMALL(2) / O402400000000, O000000000000 /
  633. C     DATA LARGE(1), LARGE(2) / O376777777777, O777777777777 /
  634. C     DATA RIGHT(1), RIGHT(2) / O604400000000, O000000000000 /
  635. C     DATA DIVER(1), DIVER(2) / O606400000000, O000000000000 /
  636. C     DATA LOG10(1), LOG10(2) / O776464202324, O117571775714 /
  637. C
  638. C     MACHINE CONSTANTS FOR THE HP 2100
  639. C     THREE WORD DOUBLE PRECISION OPTION WITH FTN4
  640. C
  641. C     DATA SMALL(1), SMALL(2), SMALL(3) / 40000B,       0,       1 /
  642. C     DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B /
  643. C     DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B,       0,    265B /
  644. C     DATA DIVER(1), DIVER(2), DIVER(3) / 40000B,       0,    276B /
  645. C     DATA LOG10(1), LOG10(2), LOG10(3) / 46420B,  46502B,  77777B /
  646. C
  647. C     MACHINE CONSTANTS FOR THE HP 2100
  648. C     FOUR WORD DOUBLE PRECISION OPTION WITH FTN4
  649. C
  650. C     DATA SMALL(1), SMALL(2) /  40000B,       0 /
  651. C     DATA SMALL(3), SMALL(4) /       0,       1 /
  652. C     DATA LARGE(1), LARGE(2) /  77777B, 177777B /
  653. C     DATA LARGE(3), LARGE(4) / 177777B, 177776B /
  654. C     DATA RIGHT(1), RIGHT(2) /  40000B,       0 /
  655. C     DATA RIGHT(3), RIGHT(4) /       0,    225B /
  656. C     DATA DIVER(1), DIVER(2) /  40000B,       0 /
  657. C     DATA DIVER(3), DIVER(4) /       0,    227B /
  658. C     DATA LOG10(1), LOG10(2) /  46420B,  46502B /
  659. C     DATA LOG10(3), LOG10(4) /  76747B, 176377B /
  660. C
  661. C     MACHINE CONSTANTS FOR THE HP 9000
  662. C
  663. C     DATA SMALL(1), SMALL(2) / 00040000000B, 00000000000B /
  664. C     DATA LARGE(1), LARGE(2) / 17737777777B, 37777777777B /
  665. C     DATA RIGHT(1), RIGHT(2) / 07454000000B, 00000000000B /
  666. C     DATA DIVER(1), DIVER(2) / 07460000000B, 00000000000B /
  667. C     DATA LOG10(1), LOG10(2) / 07764642023B, 12047674777B /
  668. C
  669. C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
  670. C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND
  671. C     THE PERKIN ELMER (INTERDATA) 7/32.
  672. C
  673. C     DATA SMALL(1), SMALL(2) / Z00100000, Z00000000 /
  674. C     DATA LARGE(1), LARGE(2) / Z7FFFFFFF, ZFFFFFFFF /
  675. C     DATA RIGHT(1), RIGHT(2) / Z33100000, Z00000000 /
  676. C     DATA DIVER(1), DIVER(2) / Z34100000, Z00000000 /
  677. C     DATA LOG10(1), LOG10(2) / Z41134413, Z509F79FF /
  678. C
  679. C     MACHINE CONSTANTS FOR THE IBM PC
  680. C     ASSUMES THAT ALL ARITHMETIC IS DONE IN DOUBLE PRECISION
  681. C     ON 8088, I.E., NOT IN 80 BIT FORM FOR THE 8087.
  682. C
  683. C     DATA SMALL(1)           / 2.23D-308              /
  684. C     DATA LARGE(1)           / 1.79D+308              /
  685. C     DATA RIGHT(1)           / 1.11D-16               /
  686. C     DATA DIVER(1)           / 2.22D-16               /
  687. C     DATA LOG10(1)           / 0.301029995663981195D0 /
  688. C
  689. C     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR)
  690. C
  691. C     DATA SMALL(1), SMALL(2) / "033400000000, "000000000000 /
  692. C     DATA LARGE(1), LARGE(2) / "377777777777, "344777777777 /
  693. C     DATA RIGHT(1), RIGHT(2) / "113400000000, "000000000000 /
  694. C     DATA DIVER(1), DIVER(2) / "114400000000, "000000000000 /
  695. C     DATA LOG10(1), LOG10(2) / "177464202324, "144117571776 /
  696. C
  697. C     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR)
  698. C
  699. C     DATA SMALL(1), SMALL(2) / "000400000000, "000000000000 /
  700. C     DATA LARGE(1), LARGE(2) / "377777777777, "377777777777 /
  701. C     DATA RIGHT(1), RIGHT(2) / "103400000000, "000000000000 /
  702. C     DATA DIVER(1), DIVER(2) / "104400000000, "000000000000 /
  703. C     DATA LOG10(1), LOG10(2) / "177464202324, "476747767461 /
  704. C
  705. C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
  706. C     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
  707. C
  708. C     DATA SMALL(1), SMALL(2) /    8388608,           0 /
  709. C     DATA LARGE(1), LARGE(2) / 2147483647,          -1 /
  710. C     DATA RIGHT(1), RIGHT(2) /  612368384,           0 /
  711. C     DATA DIVER(1), DIVER(2) /  620756992,           0 /
  712. C     DATA LOG10(1), LOG10(2) / 1067065498, -2063872008 /
  713. C
  714. C     DATA SMALL(1), SMALL(2) / O00040000000, O00000000000 /
  715. C     DATA LARGE(1), LARGE(2) / O17777777777, O37777777777 /
  716. C     DATA RIGHT(1), RIGHT(2) / O04440000000, O00000000000 /
  717. C     DATA DIVER(1), DIVER(2) / O04500000000, O00000000000 /
  718. C     DATA LOG10(1), LOG10(2) / O07746420232, O20476747770 /
  719. C
  720. C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
  721. C     16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
  722. C
  723. C     DATA SMALL(1), SMALL(2) /    128,      0 /
  724. C     DATA SMALL(3), SMALL(4) /      0,      0 /
  725. C     DATA LARGE(1), LARGE(2) /  32767,     -1 /
  726. C     DATA LARGE(3), LARGE(4) /     -1,     -1 /
  727. C     DATA RIGHT(1), RIGHT(2) /   9344,      0 /
  728. C     DATA RIGHT(3), RIGHT(4) /      0,      0 /
  729. C     DATA DIVER(1), DIVER(2) /   9472,      0 /
  730. C     DATA DIVER(3), DIVER(4) /      0,      0 /
  731. C     DATA LOG10(1), LOG10(2) /  16282,   8346 /
  732. C     DATA LOG10(3), LOG10(4) / -31493, -12296 /
  733. C
  734. C     DATA SMALL(1), SMALL(2) / O000200, O000000 /
  735. C     DATA SMALL(3), SMALL(4) / O000000, O000000 /
  736. C     DATA LARGE(1), LARGE(2) / O077777, O177777 /
  737. C     DATA LARGE(3), LARGE(4) / O177777, O177777 /
  738. C     DATA RIGHT(1), RIGHT(2) / O022200, O000000 /
  739. C     DATA RIGHT(3), RIGHT(4) / O000000, O000000 /
  740. C     DATA DIVER(1), DIVER(2) / O022400, O000000 /
  741. C     DATA DIVER(3), DIVER(4) / O000000, O000000 /
  742. C     DATA LOG10(1), LOG10(2) / O037632, O020232 /
  743. C     DATA LOG10(3), LOG10(4) / O102373, O147770 /
  744. C
  745. C     MACHINE CONSTANTS FOR THE SUN
  746. C
  747. C     DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' /
  748. C     DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' /
  749. C     DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' /
  750. C     DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' /
  751. C     DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' /
  752. C
  753. C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES FTN COMPILER
  754. C
  755. C     DATA SMALL(1), SMALL(2) / O000040000000, O000000000000 /
  756. C     DATA LARGE(1), LARGE(2) / O377777777777, O777777777777 /
  757. C     DATA RIGHT(1), RIGHT(2) / O170540000000, O000000000000 /
  758. C     DATA DIVER(1), DIVER(2) / O170640000000, O000000000000 /
  759. C     DATA LOG10(1), LOG10(2) / O177746420232, O411757177572 /
  760. C
  761. C     MACHINE CONSTANTS FOR VAX 11/780
  762. C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
  763. C     THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS
  764. C     THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS
  765. C
  766. C     DATA SMALL(1), SMALL(2) /        128,           0 /
  767. C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /
  768. C     DATA RIGHT(1), RIGHT(2) /       9344,           0 /
  769. C     DATA DIVER(1), DIVER(2) /       9472,           0 /
  770. C     DATA LOG10(1), LOG10(2) /  546979738,  -805796613 /
  771. C
  772. C     DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 /
  773. C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
  774. C     DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 /
  775. C     DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 /
  776. C     DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB /
  777. C
  778. C     MACHINE CONSTANTS FOR VAX 11/780 (G-FLOATING)
  779. C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
  780. C     THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS
  781. C     THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS
  782. C
  783. C     DATA SMALL(1), SMALL(2) /         16,           0 /
  784. C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /
  785. C     DATA RIGHT(1), RIGHT(2) /      15552,           0 /
  786. C     DATA DIVER(1), DIVER(2) /      15568,           0 /
  787. C     DATA LOG10(1), LOG10(2) /  1142112243, 2046775455 /
  788. C
  789. C     DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 /
  790. C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
  791. C     DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 /
  792. C     DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 /
  793. C     DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F /
  794. C
  795. C
  796. C***FIRST EXECUTABLE STATEMENT  D1MACH
  797. C
  798. C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
  799. C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
  800. C        NOTE ADDED BY F. ROMANI  7/11/89
  801. C
  802. C     IF (I .LT. 1  .OR.  I .GT. 5)
  803. C    1   CALL XERROR ('D1MACH -- I OUT OF BOUNDS', 25, 1, 2)
  804. C
  805.       D1MACH = DMACH(I)
  806.       RETURN
  807. C
  808.       END
  809. *DECK DQXGS
  810.       SUBROUTINE DQXGS(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,IER,
  811.      1   LIMIT,LENIW,LENW,LAST,IWORK,WORK)
  812. C
  813. C            THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
  814. C            DEFINITE INTEGRAL  I = INTEGRAL OF F OVER (A,B),
  815. C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
  816. C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
  817. C
  818. C        PARAMETERS
  819. C         ON ENTRY
  820. C            F      - DOUBLE PRECISION
  821. C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  822. C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  823. C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  824. C
  825. C            A      - DOUBLE PRECISION
  826. C                     LOWER LIMIT OF INTEGRATION
  827. C
  828. C            B      - DOUBLE PRECISION
  829. C                     UPPER LIMIT OF INTEGRATION
  830. C
  831. C            EPSABS - DOUBLE PRECISION
  832. C                     ABSOLUTE ACCURACY REQUESTED
  833. C            EPSREL - DOUBLE PRECISION
  834. C                     RELATIVE ACCURACY REQUESTED
  835. C                     IF  EPSABS.LE.0
  836. C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
  837. C                     THE ROUTINE WILL END WITH IER = 6.
  838. C
  839. C         ON RETURN
  840. C            RESULT - DOUBLE PRECISION
  841. C                     APPROXIMATION TO THE INTEGRAL
  842. C
  843. C            ABSERR - DOUBLE PRECISION
  844. C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
  845. C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
  846. C
  847. C            IER    - INTEGER
  848. C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
  849. C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
  850. C                             ACCURACY HAS BEEN ACHIEVED.
  851. C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE
  852. C                             THE ESTIMATES FOR INTEGRAL AND ERROR ARE
  853. C                             LESS RELIABLE. IT IS ASSUMED THAT THE
  854. C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
  855. C            ERROR MESSAGES
  856. C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
  857. C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB-
  858. C                             DIVISIONS BY INCREASING THE VALUE OF LIMIT
  859. C                             (AND TAKING THE ACCORDING DIMENSION
  860. C                             ADJUSTMENTS INTO ACCOUNT. HOWEVER, IF
  861. C                             THIS YIELDS NO IMPROVEMENT IT IS ADVISED
  862. C                             TO ANALYZE THE INTEGRAND IN ORDER TO
  863. C                             DETERMINE THE INTEGRATION DIFFICULTIES. IF
  864. C                             THE POSITION OF A LOCAL DIFFICULTY CAN BE
  865. C                             DETERMINED (E.G. SINGULARITY,
  866. C                             DISCONTINUITY WITHIN THE INTERVAL) ONE
  867. C                             WILL PROBABLY GAIN FROM SPLITTING UP THE
  868. C                             INTERVAL AT THIS POINT AND CALLING THE
  869. C                             INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
  870. C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
  871. C                             SHOULD BE USED, WHICH IS DESIGNED FOR
  872. C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.
  873. C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS DETEC-
  874. C                             TED, WHICH PREVENTS THE REQUESTED
  875. C                             TOLERANCE FROM BEING ACHIEVED.
  876. C                             THE ERROR MAY BE UNDER-ESTIMATED.
  877. C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR
  878. C                             OCCURS AT SOME POINTS OF THE INTEGRATION
  879. C                             INTERVAL.
  880. C                         = 4 THE ALGORITHM DOES NOT CONVERGE.
  881. C                             ROUNDOFF ERROR IS DETECTED IN THE
  882. C                             EXTRAPOLATION TABLE. IT IS PRESUMED THAT
  883. C                             THE REQUESTED TOLERANCE CANNOT BE
  884. C                             ACHIEVED, AND THAT THE RETURNED RESULT IS
  885. C                             THE BEST WHICH CAN BE OBTAINED.
  886. C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
  887. C                             SLOWLY CONVERGENT. IT MUST BE NOTED THAT
  888. C                             DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
  889. C                             OF IER.
  890. C                         = 6 THE INPUT IS INVALID, BECAUSE
  891. C                             (EPSABS.LE.0 AND
  892. C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28)
  893. C                             OR LIMIT.LT.1 OR LENW.LT.LIMIT*46 OR
  894. C                                LENIW .LT. LIMIT*3
  895. C                             RESULT, ABSERR, LAST ARE SET TO
  896. C                             ZERO.EXCEPT WHEN LIMIT OR LENW OR LENIW
  897. C                             IS INVALID, IWORK(1), WORK(LIMIT*2+1) AND
  898. C                             WORK(LIMIT*3+1) ARE SET TO ZERO, WORK(1)
  899. C                             IS SET TO A AND WORK(LIMIT+1) TO B.
  900. C
  901. C         DIMENSIONING PARAMETERS
  902. C            LIMIT - INTEGER
  903. C                    LIMIT DETERMINES THE MAXIMUM NUMBER OF SUBINTERVALS
  904. C                    IN THE PARTITION OF THE GIVEN INTEGRATION INTERVAL
  905. C                    (A,B), LIMIT.GE.1.
  906. C                    IF LIMIT.LT.1, THE ROUTINE WILL END WITH IER = 6.
  907. C
  908. C            LENW  - INTEGER
  909. C                    DIMENSIONING PARAMETER FOR WORK
  910. C                    LENW MUST BE AT LEAST LIMIT*46
  911. C                    IF LENW.LT.LIMIT*46, THE ROUTINE WILL END
  912. C                    WITH IER = 6.
  913. C
  914. C            LENIW - INTEGER
  915. C                    DIMENSIONING PARAMETER FOR IWORK
  916. C                    LENIW MUST BE AT LEAST LIMIT*3
  917. C                    IF LENW.LT.LIMIT*3, THE ROUTINE WILL END
  918. C                    WITH IER = 6.
  919. C
  920. C            LAST  - INTEGER
  921. C                    ON RETURN, LAST EQUALS THE NUMBER OF SUBINTERVALS
  922. C                    PRODUCED IN THE SUBDIVISION PROCESS, DETERMINES THE
  923. C                    NUMBER OF SIGNIFICANT ELEMENTS ACTUALLY IN THE WORK
  924. C                    ARRAYS.
  925. C
  926. C         WORK ARRAYS
  927. C            IWORK - INTEGER
  928. C                    VECTOR OF DIMENSION AT LEAST 3*LIMIT, THE FIRST K
  929. C                    ELEMENTS OF WHICH CONTAIN POINTERS
  930. C                    TO THE ERROR ESTIMATES OVER THE SUBINTERVALS
  931. C                    SUCH THAT WORK(LIMIT*3+IWORK(1)),... ,
  932. C                    WORK(LIMIT*3+IWORK(K)) FORM A DECREASING
  933. C                    SEQUENCE, WITH K = LAST IF LAST.LE.(LIMIT/2+2),
  934. C                    AND K = LIMIT+1-LAST OTHERWISE
  935. C
  936. C            WORK  - DOUBLE PRECISION
  937. C                    VECTOR OF DIMENSION AT LEAST LENW
  938. C                    ON RETURN
  939. C                    WORK(1), ..., WORK(LAST) CONTAIN THE LEFT
  940. C                     END-POINTS OF THE SUBINTERVALS IN THE
  941. C                     PARTITION OF (A,B),
  942. C                    WORK(LIMIT+1), ..., WORK(LIMIT+LAST) CONTAIN
  943. C                     THE RIGHT END-POINTS,
  944. C                    WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) CONTAIN
  945. C                     THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS,
  946. C                    WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST)
  947. C                     CONTAIN THE ERROR ESTIMATES.
  948. C                    WORK(LIMIT*4+1), ... IS THE AREA RESERVED TO STORE
  949. C                     FUNCTIONAL VALUES.
  950. C***ROUTINES CALLED  DQXGSE,XERROR
  951. C
  952. C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
  953. C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
  954. C
  955.       DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK
  956.       INTEGER IER,IWORK,LAST,LENIW,LENW,LIMIT,LVL,L1,L2,L3,L4,L5
  957. C
  958.       DIMENSION IWORK(LENIW),WORK(LENW)
  959. C
  960.       EXTERNAL F
  961. C
  962. C         CHECK VALIDITY OF LIMIT,LENIW AND LENW.
  963. C
  964. C***FIRST EXECUTABLE STATEMENT  DQXGS
  965.       IER = 6
  966.       LAST = 0
  967.       RESULT = 0.0D+00
  968.       ABSERR = 0.0D+00
  969.       IF(LIMIT.LT.1.OR.LENIW.LT.LIMIT*3.OR.LENW.LT.LIMIT*46) GO TO 10
  970. C
  971. C         PREPARE CALL FOR DQXGSE.
  972. C
  973.       L1 = LIMIT+1
  974.       L2 = LIMIT+L1
  975.       L3 = LIMIT+L2
  976.       L4 = LIMIT+L3
  977.       L5 = 21*LIMIT+L3
  978. C
  979.       CALL DQXGSE(F,A,B,EPSABS,EPSREL,LIMIT,RESULT,ABSERR,
  980.      *  IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST,
  981.      *  WORK(L4),WORK(L5),IWORK(L1),IWORK(L2))
  982. C
  983. C         CALL ERROR HANDLER IF NECESSARY.
  984. C
  985.       LVL = 0
  986. 10    IF(IER.EQ.6) LVL = 1
  987. C
  988. C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
  989. C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
  990. C
  991. C     IF(IER.NE.0)CALL XERROR('ABNORMAL RETURN FROM DQXGS  ',28,IER,LVL)
  992.       RETURN
  993.       END
  994. *DECK DQXGSE
  995.       SUBROUTINE DQXGSE(F,A,B,EPSABS,EPSREL,LIMIT,RESULT,ABSERR,
  996.      *   IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST,VALP,VALN,LP,LN)
  997. C
  998. C            THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
  999. C            DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B),
  1000. C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
  1001. C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
  1002. C
  1003. C        PARAMETERS
  1004. C         ON ENTRY
  1005. C            F      - DOUBLE PRECISION
  1006. C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  1007. C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  1008. C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  1009. C
  1010. C            A      - DOUBLE PRECISION
  1011. C                     LOWER LIMIT OF INTEGRATION
  1012. C
  1013. C            B      - DOUBLE PRECISION
  1014. C                     UPPER LIMIT OF INTEGRATION
  1015. C
  1016. C            EPSABS - DOUBLE PRECISION
  1017. C                     ABSOLUTE ACCURACY REQUESTED
  1018. C            EPSREL - DOUBLE PRECISION
  1019. C                     RELATIVE ACCURACY REQUESTED
  1020. C                     IF  EPSABS.LE.0
  1021. C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
  1022. C                     THE ROUTINE WILL END WITH IER = 6.
  1023. C
  1024. C            LIMIT  - INTEGER
  1025. C                     GIVES AN UPPERBOUND ON THE NUMBER OF SUBINTERVALS
  1026. C                     IN THE PARTITION OF (A,B)
  1027. C
  1028. C         ON RETURN
  1029. C            RESULT - DOUBLE PRECISION
  1030. C                     APPROXIMATION TO THE INTEGRAL
  1031. C
  1032. C            ABSERR - DOUBLE PRECISION
  1033. C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
  1034. C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
  1035. C
  1036. C            IER    - INTEGER
  1037. C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
  1038. C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
  1039. C                             ACCURACY HAS BEEN ACHIEVED.
  1040. C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE
  1041. C                             THE ESTIMATES FOR INTEGRAL AND ERROR ARE
  1042. C                             LESS RELIABLE. IT IS ASSUMED THAT THE
  1043. C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
  1044. C            ERROR MESSAGES
  1045. C                         = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
  1046. C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB-
  1047. C                             DIVISIONS BY INCREASING THE VALUE OF LIMIT
  1048. C                             (AND TAKING THE ACCORDING DIMENSION
  1049. C                             ADJUSTMENTS INTO ACCOUNT). HOWEVER, IF
  1050. C                             THIS YIELDS NO IMPROVEMENT IT IS ADVISED
  1051. C                             TO ANALYZE THE INTEGRAND IN ORDER TO
  1052. C                             DETERMINE THE INTEGRATION DIFFICULTIES. IF
  1053. C                             THE POSITION OF A LOCAL DIFFICULTY CAN BE
  1054. C                             DETERMINED (E.G. SINGULARITY,
  1055. C                             DISCONTINUITY WITHIN THE INTERVAL) ONE
  1056. C                             WILL PROBABLY GAIN FROM SPLITTING UP THE
  1057. C                             INTERVAL AT THIS POINT AND CALLING THE
  1058. C                             INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
  1059. C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
  1060. C                             SHOULD BE USED, WHICH IS DESIGNED FOR
  1061. C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.
  1062. C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS DETEC-
  1063. C                             TED, WHICH PREVENTS THE REQUESTED
  1064. C                             TOLERANCE FROM BEING ACHIEVED.
  1065. C                             THE ERROR MAY BE UNDER-ESTIMATED.
  1066. C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR
  1067. C                             OCCURS AT SOME POINTS OF THE INTEGRATION
  1068. C                             INTERVAL.
  1069. C                         = 4 THE ALGORITHM DOES NOT CONVERGE.
  1070. C                             ROUNDOFF ERROR IS DETECTED IN THE
  1071. C                             EXTRAPOLATION TABLE.
  1072. C                             IT IS PRESUMED THAT THE REQUESTED
  1073. C                             TOLERANCE CANNOT BE ACHIEVED, AND THAT THE
  1074. C                             RETURNED RESULT IS THE BEST WHICH CAN BE
  1075. C                             OBTAINED.
  1076. C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
  1077. C                             SLOWLY CONVERGENT. IT MUST BE NOTED THAT
  1078. C                             DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
  1079. C                             OF IER.
  1080. C                         = 6 THE INPUT IS INVALID, BECAUSE
  1081. C                             EPSABS.LE.0 AND
  1082. C                             EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28).
  1083. C                             RESULT, ABSERR, LAST, RLIST(1),
  1084. C                             IORD(1) AND ELIST(1) ARE SET TO ZERO.
  1085. C                             ALIST(1) AND BLIST(1) ARE SET TO A AND B
  1086. C                             RESPECTIVELY.
  1087. C
  1088. C            ALIST  - DOUBLE PRECISION
  1089. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  1090. C                      LAST  ELEMENTS OF WHICH ARE THE LEFT END POINTS
  1091. C                     OF THE SUBINTERVALS IN THE PARTITION OF THE
  1092. C                     GIVEN INTEGRATION RANGE (A,B)
  1093. C
  1094. C            BLIST  - DOUBLE PRECISION
  1095. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  1096. C                      LAST  ELEMENTS OF WHICH ARE THE RIGHT END POINTS
  1097. C                     OF THE SUBINTERVALS IN THE PARTITION OF THE GIVEN
  1098. C                     INTEGRATION RANGE (A,B)
  1099. C
  1100. C            RLIST  - DOUBLE PRECISION
  1101. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  1102. C                      LAST  ELEMENTS OF WHICH ARE THE INTEGRAL
  1103. C                     APPROXIMATIONS ON THE SUBINTERVALS
  1104. C
  1105. C            ELIST  - DOUBLE PRECISION
  1106. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  1107. C                      LAST  ELEMENTS OF WHICH ARE THE MODULI OF THE
  1108. C                     ABSOLUTE ERROR ESTIMATES ON THE SUBINTERVALS
  1109. C
  1110. C            IORD   - INTEGER
  1111. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST K
  1112. C                     ELEMENTS OF WHICH ARE POINTERS TO THE
  1113. C                     ERROR ESTIMATES OVER THE SUBINTERVALS,
  1114. C                     SUCH THAT ELIST(IORD(1)), ..., ELIST(IORD(K))
  1115. C                     FORM A DECREASING SEQUENCE, WITH K = LAST
  1116. C                     IF LAST.LE.(LIMIT/2+2), AND K = LIMIT+1-LAST
  1117. C                     OTHERWISE
  1118. C
  1119. C            LAST   - INTEGER
  1120. C                     NUMBER OF SUBINTERVALS ACTUALLY PRODUCED IN THE
  1121. C                     SUBDIVISION PROCESS
  1122. C
  1123. C            VALP   - DOUBLE PRECISION
  1124. C            VALN     ARRAYS OF DIMENSION AT LEAST (21,LIMIT) USED TO
  1125. C                     SAVE THE FUNCTIONAL VALUES
  1126. C
  1127. C            LP     - INTEGER
  1128. C            LN       VECTORS OF DIMENSION AT LEAST LIMIT, USED TO
  1129. C                     STORE THE ACTUAL NUMBER OF FUNCTIONAL VALUES
  1130. C                     SAVED IN THE CORRESPONDING COLUMN
  1131. C                     OF VALP,VALN
  1132. C
  1133. C***ROUTINES CALLED  F,D1MACH,DQELG,DQXLQM,DQPSRT,DQXRRD,DQXCPY
  1134. C
  1135.       DOUBLE PRECISION A,ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,
  1136.      *  A2,B,BLIST,B1,B2,CORREC,DABS,DEFABS,DEFAB1,DEFAB2,D1MACH,DMAX1,
  1137.      *  DRES,ELIST,EPMACH,EPSABS,EPSREL,ERLARG,ERLAST,ERRBND,ERRMAX,
  1138.      *  ERROR1,ERROR2,ERRO12,ERRSUM,ERTEST,F,OFLOW,RESABS,RESEPS,RESULT,
  1139.      *  RES3LA,RLIST,RLIST2,SMALL,UFLOW,
  1140.      *  VALP,VALN,VP1,VP2,VN1,VN2
  1141.       INTEGER ID,IER,IERRO,IORD,IROFF1,IROFF2,IROFF3,JUPBND,K,KSGN,
  1142.      *  KTMIN,LAST,LIMIT,MAXERR,NRES,NRMAX,NUMRL2,
  1143.      *  LP,LN,LP1,LP2,LN1,LN2
  1144.       LOGICAL EXTRAP,NOEXT
  1145.       DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT),
  1146.      * RES3LA(3),RLIST(LIMIT),RLIST2(52),
  1147.      * VALP(21,LIMIT),VALN(21,LIMIT),LP(LIMIT),LN(LIMIT),
  1148.      * VP1(21),VP2(21),VN1(21),VN2(21)
  1149.       EXTERNAL F
  1150. C
  1151. C            MACHINE DEPENDENT CONSTANTS
  1152. C            ---------------------------
  1153. C
  1154. C          EPMACH IS THE LARGEST RELATIVE SPACING.
  1155. C          UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
  1156. C          OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
  1157. C
  1158. C***FIRST EXECUTABLE STATEMENT  DQXGSE
  1159.       EPMACH = D1MACH(4)
  1160. C
  1161. C            TEST ON VALIDITY OF PARAMETERS
  1162. C            ------------------------------
  1163.       IER = 0
  1164.       LAST = 0
  1165.       RESULT = 0.0D+00
  1166.       ABSERR = 0.0D+00
  1167.       ALIST(1) = A
  1168.       BLIST(1) = B
  1169.       RLIST(1) = 0.0D+00
  1170.       ELIST(1) = 0.0D+00
  1171.       IF(EPSABS.LE.0.0D+00.AND.EPSREL.LT.DMAX1(0.5D+02*EPMACH,0.5D-28))
  1172.      *   IER = 6
  1173.       IF(IER.EQ.6) GO TO 999
  1174. C
  1175. C           FIRST APPROXIMATION TO THE INTEGRAL
  1176. C           -----------------------------------
  1177. C
  1178.       UFLOW = D1MACH(1)
  1179.       OFLOW = D1MACH(2)
  1180.       IERRO = 0
  1181.       LP(1)=1
  1182.       LN(1)=1
  1183.       VALP(1,1)=F((A+B)*0.5D0)
  1184.       VALN(1,1)=VALP(1,1)
  1185.       CALL DQXLQM(F,A,B,RESULT,ABSERR,DEFABS,RESABS,
  1186.      *           VALP(1,1),VALN(1,1),LP(1),LN(1), 2  )
  1187. C
  1188. C           TEST ON ACCURACY.
  1189. C
  1190.       DRES = DABS(RESULT)
  1191.       ERRBND = DMAX1(EPSABS,EPSREL*DRES)
  1192.       LAST = 1
  1193.       RLIST(1) = RESULT
  1194.       ELIST(1) = ABSERR
  1195.       IORD(1) = 1
  1196.       IF(ABSERR.LE.1.0D+02*EPMACH*DEFABS.AND.ABSERR.GT.ERRBND) IER = 2
  1197.       IF(LIMIT.EQ.1) IER = 1
  1198.       IF(IER.NE.0.OR.(ABSERR.LE.ERRBND.AND.ABSERR.NE.RESABS).OR.
  1199.      *  ABSERR.EQ.0.0D+00) GO TO 999
  1200. C
  1201. C           INITIALIZATION
  1202. C           --------------
  1203. C
  1204.       RLIST2(1) = RESULT
  1205.       ERRMAX = ABSERR
  1206.       MAXERR = 1
  1207.       AREA = RESULT
  1208.       ERRSUM = ABSERR
  1209.       ABSERR = OFLOW
  1210.       NRMAX = 1
  1211.       NRES = 0
  1212.       NUMRL2 = 2
  1213.       KTMIN = 0
  1214.       EXTRAP = .FALSE.
  1215.       NOEXT = .FALSE.
  1216.       IROFF1 = 0
  1217.       IROFF2 = 0
  1218.       IROFF3 = 0
  1219.       KSGN = -1
  1220.       IF(DRES.GE.(0.1D+01-0.5D+02*EPMACH)*DEFABS) KSGN = 1
  1221. C
  1222. C           MAIN DO-LOOP
  1223. C           ------------
  1224. C
  1225.       DO 90 LAST = 2,LIMIT
  1226. C
  1227. C           BISECT THE SUBINTERVAL WITH THE NRMAX-TH LARGEST ERROR
  1228. C           ESTIMATE.
  1229. C
  1230.         A1 = ALIST(MAXERR)
  1231.         B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR))
  1232.         A2 = B1
  1233.         B2 = BLIST(MAXERR)
  1234.         ERLAST = ERRMAX
  1235.         CALL DQXRRD(F,VALN(1,MAXERR),LN(MAXERR),B1,A1,VN1,VP1,LN1,LP1)
  1236.         CALL DQXLQM(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1,VP1,VN1,LP1,LN1,
  1237.      *              2)
  1238.         CALL DQXRRD(F,VALP(1,MAXERR),LP(MAXERR),A2,B2,VP2,VN2,LP2,LN2)
  1239.         CALL DQXLQM(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2,VP2,VN2,LP2,LN2,
  1240.      *              2)
  1241. C
  1242. C           IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
  1243. C           AND ERROR AND TEST FOR ACCURACY.
  1244. C
  1245.         AREA12 = AREA1+AREA2
  1246.         ERRO12 = ERROR1+ERROR2
  1247.         ERRSUM = ERRSUM+ERRO12-ERRMAX
  1248.         AREA = AREA+AREA12-RLIST(MAXERR)
  1249.         IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2) GO TO 15
  1250.         IF(DABS(RLIST(MAXERR)-AREA12).GT.0.1D-04*DABS(AREA12)
  1251.      *  .OR.ERRO12.LT.0.99D+00*ERRMAX) GO TO 10
  1252.         IF(EXTRAP) IROFF2 = IROFF2+1
  1253.         IF(.NOT.EXTRAP) IROFF1 = IROFF1+1
  1254.    10   IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF3 = IROFF3+1
  1255.    15   RLIST(MAXERR) = AREA1
  1256.         RLIST(LAST) = AREA2
  1257.         ERRBND = DMAX1(EPSABS,EPSREL*DABS(AREA))
  1258. C
  1259. C           TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
  1260. C
  1261.         IF(IROFF1+IROFF2.GE.10.OR.IROFF3.GE.20) IER = 2
  1262.         IF(IROFF2.GE.5) IERRO = 3
  1263. C
  1264. C           SET ERROR FLAG IN THE CASE THAT THE NUMBER OF SUBINTERVALS
  1265. C           EQUALS LIMIT.
  1266. C
  1267.         IF(LAST.EQ.LIMIT) IER = 1
  1268. C
  1269. C           SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
  1270. C           AT A POINT OF THE INTEGRATION RANGE.
  1271. C
  1272.         IF(DMAX1(DABS(A1),DABS(B2)).LE.(0.1D+01+0.1D+03*EPMACH)*
  1273.      *  (DABS(A2)+0.1D+04*UFLOW)) IER = 4
  1274. C
  1275. C           APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
  1276. C
  1277.         IF(ERROR2.GT.ERROR1) GO TO 20
  1278.         ALIST(LAST) = A2
  1279.         BLIST(MAXERR) = B1
  1280.         BLIST(LAST) = B2
  1281.         ELIST(MAXERR) = ERROR1
  1282.         ELIST(LAST) = ERROR2
  1283.         CALL DQXCPY(VALP(1,MAXERR),VP1,LP1)
  1284.         LP(MAXERR)=LP1
  1285.         CALL DQXCPY(VALN(1,MAXERR),VN1,LN1)
  1286.         LN(MAXERR)=LN1
  1287.         CALL DQXCPY(VALP(1,LAST),VP2,LP2)
  1288.         LP(LAST)=LP2
  1289.         CALL DQXCPY(VALN(1,LAST),VN2,LN2)
  1290.         LN(LAST)=LN2
  1291.         GO TO 30
  1292.    20   ALIST(MAXERR) = A2
  1293.         ALIST(LAST) = A1
  1294.         BLIST(LAST) = B1
  1295.         RLIST(MAXERR) = AREA2
  1296.         RLIST(LAST) = AREA1
  1297.         ELIST(MAXERR) = ERROR2
  1298.         ELIST(LAST) = ERROR1
  1299.         CALL DQXCPY(VALP(1,MAXERR),VP2,LP2)
  1300.         LP(MAXERR)=LP2
  1301.         CALL DQXCPY(VALN(1,MAXERR),VN2,LN2)
  1302.         LN(MAXERR)=LN2
  1303.         CALL DQXCPY(VALP(1,LAST),VP1,LP1)
  1304.         LP(LAST)=LP1
  1305.         CALL DQXCPY(VALN(1,LAST),VN1,LN1)
  1306.         LN(LAST)=LN1
  1307. C
  1308. C           CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING
  1309. C           IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
  1310. C           WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
  1311. C
  1312.    30   CALL DQPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX)
  1313. C ***JUMP OUT OF DO-LOOP
  1314.         IF(ERRSUM.LE.ERRBND) GO TO 115
  1315. C ***JUMP OUT OF DO-LOOP
  1316.         IF(IER.NE.0) GO TO 100
  1317.         IF(LAST.EQ.2) GO TO 80
  1318.         IF(NOEXT) GO TO 90
  1319.         ERLARG = ERLARG-ERLAST
  1320.         IF(DABS(B1-A1).GT.SMALL) ERLARG = ERLARG+ERRO12
  1321.         IF(EXTRAP) GO TO 40
  1322. C
  1323. C           TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE
  1324. C           SMALLEST INTERVAL.
  1325. C
  1326.         IF(DABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90
  1327.         EXTRAP = .TRUE.
  1328.         NRMAX = 2
  1329. C
  1330. C           THE BOUND 0.3*ERTEST HAS BEEN INTRODUCED TO PERFORM A
  1331. C           MORE CAUTIOUS EXTRAPOLATION THAN IN THE ORIGINAL DQAGSE
  1332. C           ROUTINE
  1333. C
  1334.    40   IF(IERRO.EQ.3.OR.ERLARG.LE.0.3D0*ERTEST) GO TO 60
  1335. C
  1336. C           THE SMALLEST INTERVAL HAS THE LARGEST ERROR.
  1337. C           BEFORE BISECTING DECREASE THE SUM OF THE ERRORS OVER THE
  1338. C           LARGER INTERVALS (ERLARG) AND PERFORM EXTRAPOLATION.
  1339. C
  1340.         ID = NRMAX
  1341.         JUPBND = LAST
  1342.         IF(LAST.GT.(2+LIMIT/2)) JUPBND = LIMIT+3-LAST
  1343.         DO 50 K = ID,JUPBND
  1344.           MAXERR = IORD(NRMAX)
  1345.           ERRMAX = ELIST(MAXERR)
  1346. C ***JUMP OUT OF DO-LOOP
  1347.           IF(DABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90
  1348.           NRMAX = NRMAX+1
  1349.    50   CONTINUE
  1350. C
  1351. C           PERFORM EXTRAPOLATION.
  1352. C
  1353.    60   NUMRL2 = NUMRL2+1
  1354.         RLIST2(NUMRL2) = AREA
  1355.         CALL DQELG(NUMRL2,RLIST2,RESEPS,ABSEPS,RES3LA,NRES)
  1356.         KTMIN = KTMIN+1
  1357.         IF(KTMIN.GT.5.AND.ABSERR.LT.0.1D-02*ERRSUM) IER = 5
  1358.         IF(ABSEPS.GE.ABSERR) GO TO 70
  1359.         KTMIN = 0
  1360.         ABSERR = ABSEPS
  1361.         RESULT = RESEPS
  1362.         CORREC = ERLARG
  1363.         ERTEST = DMAX1(EPSABS,EPSREL*DABS(RESEPS))
  1364. C ***JUMP OUT OF DO-LOOP
  1365.         IF(ABSERR.LE.ERTEST) GO TO 100
  1366. C
  1367. C           PREPARE BISECTION OF THE SMALLEST INTERVAL.
  1368. C
  1369.    70   IF(NUMRL2.EQ.1) NOEXT = .TRUE.
  1370.         IF(IER.EQ.5) GO TO 100
  1371.         MAXERR = IORD(1)
  1372.         ERRMAX = ELIST(MAXERR)
  1373.         NRMAX = 1
  1374.         EXTRAP = .FALSE.
  1375.         SMALL = SMALL*0.5D+00
  1376.         ERLARG = ERRSUM
  1377.         GO TO 90
  1378.    80   SMALL = DABS(B-A)*0.375D+00
  1379.         ERLARG = ERRSUM
  1380.         ERTEST = ERRBND
  1381.         RLIST2(2) = AREA
  1382.    90 CONTINUE
  1383. C
  1384. C           SET FINAL RESULT AND ERROR ESTIMATE.
  1385. C           ------------------------------------
  1386. C
  1387.   100 IF(ABSERR.EQ.OFLOW) GO TO 115
  1388.       IF(IER+IERRO.EQ.0) GO TO 110
  1389.       IF(IERRO.EQ.3) ABSERR = ABSERR+CORREC
  1390.       IF(IER.EQ.0) IER = 3
  1391.       IF(RESULT.NE.0.0D+00.AND.AREA.NE.0.0D+00) GO TO 105
  1392.       IF(ABSERR.GT.ERRSUM) GO TO 115
  1393.       IF(AREA.EQ.0.0D+00) GO TO 130
  1394.       GO TO 110
  1395.   105 IF(ABSERR/DABS(RESULT).GT.ERRSUM/DABS(AREA)) GO TO 115
  1396. C
  1397. C           TEST ON DIVERGENCE.
  1398. C
  1399.   110 IF(KSGN.EQ.(-1).AND.DMAX1(DABS(RESULT),DABS(AREA)).LE.
  1400.      * DEFABS*0.1D-01) GO TO 130
  1401.       IF(0.1D-01.GT.(RESULT/AREA).OR.(RESULT/AREA).GT.0.1D+03
  1402.      * .OR.ERRSUM.GT.DABS(AREA)) IER = 6
  1403.       GO TO 130
  1404. C
  1405. C           COMPUTE GLOBAL INTEGRAL SUM.
  1406. C
  1407.   115 RESULT = 0.0D+00
  1408.       DO 120 K = 1,LAST
  1409.          RESULT = RESULT+RLIST(K)
  1410.   120 CONTINUE
  1411.       ABSERR = ERRSUM
  1412.   130 IF(IER.GT.2) IER = IER-1
  1413.   999 RETURN
  1414.       END
  1415. *DECK DQXG
  1416.       SUBROUTINE DQXG(F,A,B,EPSABS,EPSREL,KEY,RESULT,ABSERR,IER,
  1417.      1   LIMIT,LENIW,LENW,LAST,IWORK,WORK)
  1418. C
  1419. C            THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
  1420. C            DEFINITE INTEGRAL  I = INTEGRAL OF F OVER (A,B),
  1421. C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
  1422. C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
  1423. C
  1424. C        PARAMETERS
  1425. C         ON ENTRY
  1426. C            F      - DOUBLE PRECISION
  1427. C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  1428. C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  1429. C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  1430. C
  1431. C            A      - DOUBLE PRECISION
  1432. C                     LOWER LIMIT OF INTEGRATION
  1433. C
  1434. C            B      - DOUBLE PRECISION
  1435. C                     UPPER LIMIT OF INTEGRATION
  1436. C
  1437. C            EPSABS - DOUBLE PRECISION
  1438. C                     ABSOLUTE ACCURACY REQUESTED
  1439. C            EPSREL - DOUBLE PRECISION
  1440. C                     RELATIVE ACCURACY REQUESTED
  1441. C                     IF  EPSABS.LE.0
  1442. C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
  1443. C                     THE ROUTINE WILL END WITH IER = 6.
  1444. C
  1445. C            KEY    - INTEGER
  1446. C                     KEY FOR CHOICE OF LOCAL INTEGRATION RULE
  1447. C                     RMS FORMULAS ARE USED WITH
  1448. C                      13 - 19               POINTS IF KEY.LT.1,
  1449. C                      13 - 19 - (27)        POINTS IF KEY = 1,
  1450. C                      13 - 19 - (27) - (41) POINTS IF KEY = 2,
  1451. C                           19 -  27  - (41) POINTS IF KEY = 3,
  1452. C                                 27  -  41  POINTS IF KEY.GT.3.
  1453. C
  1454. C                         (RULES) USED IF THE FUNCTION APPEARS
  1455. C                         ENOUGH REGULAR IN THE SUBINTERVAL
  1456. C
  1457. C         ON RETURN
  1458. C            RESULT - DOUBLE PRECISION
  1459. C                     APPROXIMATION TO THE INTEGRAL
  1460. C
  1461. C            ABSERR - DOUBLE PRECISION
  1462. C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
  1463. C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
  1464. C
  1465. C            IER    - INTEGER
  1466. C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
  1467. C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
  1468. C                             ACCURACY HAS BEEN ACHIEVED.
  1469. C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE
  1470. C                             THE ESTIMATES FOR RESULT AND ERROR ARE
  1471. C                             LESS RELIABLE. IT IS ASSUMED THAT THE
  1472. C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
  1473. C            ERROR MESSAGES
  1474. C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
  1475. C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB-
  1476. C                             DIVISIONS BY INCREASING THE VALUE OF LIMIT
  1477. C                             (AND TAKING THE ACCORDING DIMENSION
  1478. C                             ADJUSTMENTS INTO ACCOUNT. HOWEVER, IF
  1479. C                             THIS YIELDS NO IMPROVEMENT IT IS ADVISED
  1480. C                             TO ANALYZE THE INTEGRAND IN ORDER TO
  1481. C                             DETERMINE THE INTEGRATION DIFFICULTIES. IF
  1482. C                             THE POSITION OF A LOCAL DIFFICULTY CAN BE
  1483. C                             DETERMINED (E.G. SINGULARITY,
  1484. C                             DISCONTINUITY WITHIN THE INTERVAL) ONE
  1485. C                             WILL PROBABLY GAIN FROM SPLITTING UP THE
  1486. C                             INTERVAL AT THIS POINT AND CALLING THE
  1487. C                             INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
  1488. C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
  1489. C                             SHOULD BE USED, WHICH IS DESIGNED FOR
  1490. C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.
  1491. C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS DETEC-
  1492. C                             TED, WHICH PREVENTS THE REQUESTED
  1493. C                             TOLERANCE FROM BEING ACHIEVED.
  1494. C                             THE ERROR MAY BE UNDER-ESTIMATED.
  1495. C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR
  1496. C                             OCCURS AT SOME POINTS OF THE INTEGRATION
  1497. C                             INTERVAL.
  1498. C                         = 4 THE ALGORITHM DOES NOT CONVERGE.
  1499. C                             ROUNDOFF ERROR IS DETECTED IN THE
  1500. C                             EXTRAPOLATION TABLE. IT IS PRESUMED THAT
  1501. C                             THE REQUESTED TOLERANCE CANNOT BE
  1502. C                             ACHIEVED, AND THAT THE RETURNED RESULT IS
  1503. C                             THE BEST WHICH CAN BE OBTAINED.
  1504. C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
  1505. C                             SLOWLY CONVERGENT. IT MUST BE NOTED THAT
  1506. C                             DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
  1507. C                             OF IER.
  1508. C                         = 6 THE INPUT IS INVALID, BECAUSE
  1509. C                             (EPSABS.LE.0 AND
  1510. C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28)
  1511. C                             OR LIMIT.LT.1 OR LENW.LT.LIMIT*46 OR
  1512. C                                LENIW .LT. LIMIT*3
  1513. C                             RESULT, ABSERR, LAST ARE SET TO
  1514. C                             ZERO.EXCEPT WHEN LIMIT OR LENW OR LENIW
  1515. C                             IS INVALID, IWORK(1), WORK(LIMIT*2+1) AND
  1516. C                             WORK(LIMIT*3+1) ARE SET TO ZERO, WORK(1)
  1517. C                             IS SET TO A AND WORK(LIMIT+1) TO B.
  1518. C
  1519. C         DIMENSIONING PARAMETERS
  1520. C            LIMIT - INTEGER
  1521. C                    LIMIT DETERMINES THE MAXIMUM NUMBER OF SUBINTERVALS
  1522. C                    IN THE PARTITION OF THE GIVEN INTEGRATION INTERVAL
  1523. C                    (A,B), LIMIT.GE.1.
  1524. C                    IF LIMIT.LT.1, THE ROUTINE WILL END WITH IER = 6.
  1525. C
  1526. C            LENW  - INTEGER
  1527. C                    DIMENSIONING PARAMETER FOR WORK
  1528. C                    LENW MUST BE AT LEAST LIMIT*46
  1529. C                    IF LENW.LT.LIMIT*46, THE ROUTINE WILL END
  1530. C                    WITH IER = 6.
  1531. C
  1532. C            LENIW - INTEGER
  1533. C                    DIMENSIONING PARAMETER FOR IWORK
  1534. C                    LENIW MUST BE AT LEAST LIMIT*3
  1535. C                    IF LENW.LT.LIMIT*3, THE ROUTINE WILL END
  1536. C                    WITH IER = 6.
  1537. C
  1538. C            LAST  - INTEGER
  1539. C                    ON RETURN, LAST EQUALS THE NUMBER OF SUBINTERVALS
  1540. C                    PRODUCED IN THE SUBDIVISION PROCESS, DETERMINES THE
  1541. C                    NUMBER OF SIGNIFICANT ELEMENTS ACTUALLY IN THE WORK
  1542. C                    ARRAYS.
  1543. C
  1544. C         WORK ARRAYS
  1545. C            IWORK - INTEGER
  1546. C                    VECTOR OF DIMENSION AT LEAST 3*LIMIT, THE FIRST K
  1547. C                    ELEMENTS OF WHICH CONTAIN POINTERS
  1548. C                    TO THE ERROR ESTIMATES OVER THE SUBINTERVALS
  1549. C                    SUCH THAT WORK(LIMIT*3+IWORK(1)),... ,
  1550. C                    WORK(LIMIT*3+IWORK(K)) FORM A DECREASING
  1551. C                    SEQUENCE, WITH K = LAST IF LAST.LE.(LIMIT/2+2),
  1552. C                    AND K = LIMIT+1-LAST OTHERWISE
  1553. C
  1554. C            WORK  - DOUBLE PRECISION
  1555. C                    VECTOR OF DIMENSION AT LEAST LENW
  1556. C                    ON RETURN
  1557. C                    WORK(1), ..., WORK(LAST) CONTAIN THE LEFT
  1558. C                     END-POINTS OF THE SUBINTERVALS IN THE
  1559. C                     PARTITION OF (A,B),
  1560. C                    WORK(LIMIT+1), ..., WORK(LIMIT+LAST) CONTAIN
  1561. C                     THE RIGHT END-POINTS,
  1562. C                    WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) CONTAIN
  1563. C                     THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS,
  1564. C                    WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST)
  1565. C                     CONTAIN THE ERROR ESTIMATES.
  1566. C                    WORK(LIMIT*4+1), ... IS THE AREA RESERVED TO STORE
  1567. C                     FUNCTIONAL VALUES.
  1568. C***ROUTINES CALLED  DQXGE,XERROR
  1569. C
  1570. C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
  1571. C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
  1572. C
  1573. C
  1574.       DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK
  1575.       INTEGER KEY,IER,IWORK,LAST,LENIW,LENW,LIMIT,LVL,L1,L2,L3,L4,L5
  1576. C
  1577.       DIMENSION IWORK(LENIW),WORK(LENW)
  1578. C
  1579.       EXTERNAL F
  1580. C
  1581. C         CHECK VALIDITY OF LIMIT,LENIW AND LENW.
  1582. C
  1583. C***FIRST EXECUTABLE STATEMENT  DQXG
  1584.       IER = 6
  1585.       LAST = 0
  1586.       RESULT = 0.0D+00
  1587.       ABSERR = 0.0D+00
  1588.       IF(LIMIT.LT.1.OR.LENIW.LT.LIMIT*3.OR.LENW.LT.LIMIT*46) GO TO 10
  1589. C
  1590. C         PREPARE CALL FOR DQXGE.
  1591. C
  1592.       L1 = LIMIT+1
  1593.       L2 = LIMIT+L1
  1594.       L3 = LIMIT+L2
  1595.       L4 = LIMIT+L3
  1596.       L5 = 21*LIMIT+L3
  1597. C
  1598.       CALL DQXGE(F,A,B,EPSABS,EPSREL,KEY,LIMIT,RESULT,ABSERR,
  1599.      *  IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST,
  1600.      *  WORK(L4),WORK(L5),IWORK(L1),IWORK(L2))
  1601. C
  1602. C         CALL ERROR HANDLER IF NECESSARY.
  1603. C
  1604.       LVL = 0
  1605. 10    IF(IER.EQ.6) LVL = 1
  1606. C
  1607. C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
  1608. C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
  1609. C
  1610. C     IF(IER.NE.0)CALL XERROR('ABNORMAL RETURN FROM DQXG   ',28,IER,LVL)
  1611.       RETURN
  1612.       END
  1613. *DECK DQXGE
  1614.       SUBROUTINE DQXGE(F,A,B,EPSABS,EPSREL,KEY,LIMIT,RESULT,ABSERR,
  1615.      *   IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST,VALP,VALN,LP,LN)
  1616. C
  1617. C
  1618. C            THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
  1619. C            DEFINITE INTEGRAL   I = INTEGRAL OF F OVER (A,B),
  1620. C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
  1621. C            ABS(I-RESLT).LE.MAX(EPSABS,EPSREL*ABS(I)).
  1622. C
  1623. C        PARAMETERS
  1624. C         ON ENTRY
  1625. C            F      - DOUBLE PRECISION
  1626. C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  1627. C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  1628. C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  1629. C
  1630. C            A      - DOUBLE PRECISION
  1631. C                     LOWER LIMIT OF INTEGRATION
  1632. C
  1633. C            B      - DOUBLE PRECISION
  1634. C                     UPPER LIMIT OF INTEGRATION
  1635. C
  1636. C            EPSABS - DOUBLE PRECISION
  1637. C                     ABSOLUTE ACCURACY REQUESTED
  1638. C
  1639. C            EPSREL - DOUBLE PRECISION
  1640. C                     RELATIVE ACCURACY REQUESTED
  1641. C                     IF  EPSABS.LE.0
  1642. C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
  1643. C                     THE ROUTINE WILL END WITH IER = 6.
  1644. C
  1645. C            KEY    - INTEGER
  1646. C                     KEY FOR CHOICE OF LOCAL INTEGRATION RULE
  1647. C                     RMS FORMULAS ARE USED WITH
  1648. C                      13 - 19               POINTS IF KEY.LT.1,
  1649. C                      13 - 19 - (27)        POINTS IF KEY = 1,
  1650. C                      13 - 19 - (27) - (41) POINTS IF KEY = 2,
  1651. C                           19 -  27  - (41) POINTS IF KEY = 3,
  1652. C                                 27  -  41  POINTS IF KEY.GT.3.
  1653. C
  1654. C                         (RULES) USED IF THE FUNCTION APPEARS
  1655. C                         ENOUGH REGULAR IN THE SUBINTERVAL
  1656. C
  1657. C            LIMIT  - INTEGER
  1658. C                     GIVES AN UPPERBOUND ON THE NUMBER OF SUBINTERVALS
  1659. C                     IN THE PARTITION OF (A,B), LIMIT.GE.1.
  1660. C
  1661. C         ON RETURN
  1662. C            RESULT - DOUBLE PRECISION
  1663. C                     APPROXIMATION TO THE INTEGRAL
  1664. C
  1665. C            ABSERR - DOUBLE PRECISION
  1666. C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
  1667. C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
  1668. C
  1669. C            IER    - INTEGER
  1670. C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
  1671. C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
  1672. C                             ACCURACY HAS BEEN ACHIEVED.
  1673. C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE
  1674. C                             THE ESTIMATES FOR RESULT AND ERROR ARE
  1675. C                             LESS RELIABLE. IT IS ASSUMED THAT THE
  1676. C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
  1677. C            ERROR MESSAGES
  1678. C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
  1679. C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE
  1680. C                             SUBDIVISIONS BY INCREASING THE VALUE
  1681. C                             OF LIMIT.
  1682. C                             HOWEVER, IF THIS YIELDS NO IMPROVEMENT IT
  1683. C                             IS RATHER ADVISED TO ANALYZE THE INTEGRAND
  1684. C                             IN ORDER TO DETERMINE THE INTEGRATION
  1685. C                             DIFFICULTIES. IF THE POSITION OF A LOCAL
  1686. C                             DIFFICULTY CAN BE DETERMINED(E.G.
  1687. C                             SINGULARITY, DISCONTINUITY WITHIN THE
  1688. C                             INTERVAL) ONE WILL PROBABLY GAIN FROM
  1689. C                             SPLITTING UP THE INTERVAL AT THIS POINT
  1690. C                             AND CALLING THE INTEGRATOR ON THE
  1691. C                             SUBRANGES. IF POSSIBLE, AN APPROPRIATE
  1692. C                             SPECIAL-PURPOSE INTEGRATOR SHOULD BE USED
  1693. C                             WHICH IS DESIGNED FOR HANDLING THE TYPE OF
  1694. C                             DIFFICULTY INVOLVED.
  1695. C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS
  1696. C                             DETECTED, WHICH PREVENTS THE REQUESTED
  1697. C                             TOLERANCE FROM BEING ACHIEVED.
  1698. C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS
  1699. C                             AT SOME POINTS OF THE INTEGRATION
  1700. C                             INTERVAL.
  1701. C                         = 6 THE INPUT IS INVALID, BECAUSE
  1702. C                             (EPSABS.LE.0 AND
  1703. C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
  1704. C                             RESULT, ABSERR, LAST, RLIST(1) ,
  1705. C                             ELIST(1) AND IORD(1) ARE SET TO ZERO.
  1706. C                             ALIST(1) AND BLIST(1) ARE SET TO A AND B
  1707. C                             RESPECTIVELY.
  1708. C
  1709. C            ALIST   - DOUBLE PRECISION
  1710. C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  1711. C                       LAST  ELEMENTS OF WHICH ARE THE LEFT
  1712. C                      END POINTS OF THE SUBINTERVALS IN THE PARTITION
  1713. C                      OF THE GIVEN INTEGRATION RANGE (A,B)
  1714. C
  1715. C            BLIST   - DOUBLE PRECISION
  1716. C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  1717. C                       LAST  ELEMENTS OF WHICH ARE THE RIGHT
  1718. C                      END POINTS OF THE SUBINTERVALS IN THE PARTITION
  1719. C                      OF THE GIVEN INTEGRATION RANGE (A,B)
  1720. C
  1721. C            RLIST   - DOUBLE PRECISION
  1722. C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  1723. C                       LAST  ELEMENTS OF WHICH ARE THE
  1724. C                      INTEGRAL APPROXIMATIONS ON THE SUBINTERVALS
  1725. C
  1726. C            ELIST   - DOUBLE PRECISION
  1727. C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  1728. C                       LAST  ELEMENTS OF WHICH ARE THE MODULI OF THE
  1729. C                      ABSOLUTE ERROR ESTIMATES ON THE SUBINTERVALS
  1730. C
  1731. C            IORD    - INTEGER
  1732. C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST K
  1733. C                      ELEMENTS OF WHICH ARE POINTERS TO THE
  1734. C                      ERROR ESTIMATES OVER THE SUBINTERVALS,
  1735. C                      SUCH THAT ELIST(IORD(1)), ...,
  1736. C                      ELIST(IORD(K)) FORM A DECREASING SEQUENCE,
  1737. C                      WITH K = LAST IF LAST.LE.(LIMIT/2+2), AND
  1738. C                      K = LIMIT+1-LAST OTHERWISE
  1739. C
  1740. C            LAST    - INTEGER
  1741. C                      NUMBER OF SUBINTERVALS ACTUALLY PRODUCED IN THE
  1742. C                      SUBDIVISION PROCESS
  1743. C
  1744. C            VALP   - DOUBLE PRECISION
  1745. C            VALN     ARRAYS OF DIMENSION AT LEAST (21,LIMIT) USED TO
  1746. C                     SAVE THE FUNCTIONAL VALUES
  1747. C
  1748. C            LP     - INTEGER
  1749. C            LN       VECTORS OF DIMENSION AT LEAST LIMIT, USED TO
  1750. C                     STORE THE ACTUAL NUMBER OF FUNCTIONAL VALUES
  1751. C                     SAVED IN THE CORRESPONDING COLUMN
  1752. C                     OF VALP,VALN
  1753. C
  1754. C***ROUTINES CALLED  F,D1MACH,DQXLQM,DQXRRD,DQPSRT,DQXCPY
  1755. C
  1756.       DOUBLE PRECISION A,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,A2,B,
  1757.      *  BLIST,B1,B2,DABS,DEFABS,DEFAB1,DEFAB2,DMAX1,D1MACH,ELIST,EPMACH,
  1758.      *  EPSABS,EPSREL,ERRBND,ERRMAX,ERROR1,ERROR2,ERRO12,ERRSUM,F,
  1759.      *  RESABS,RESULT,RLIST,UFLOW,VALP,VALN,VP1,VP2,VN1,VN2
  1760.       INTEGER KEY,IER,IORD,IROFF1,IROFF2,K,LAST,LIMIT,MAXERR,
  1761.      *  NRMAX,LP,LN,LP1,LP2,LN1,LN2
  1762. C
  1763.       DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT),
  1764.      *  RLIST(LIMIT),VALP(21,LIMIT),VALN(21,LIMIT),LP(LIMIT),LN(LIMIT),
  1765.      * VP1(21),VP2(21),VN1(21),VN2(21)
  1766. C
  1767.       EXTERNAL F
  1768. C
  1769. C            MACHINE DEPENDENT CONSTANTS
  1770. C            ---------------------------
  1771. C
  1772. C          EPMACH IS THE LARGEST RELATIVE SPACING.
  1773. C          UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
  1774. C
  1775. C***FIRST EXECUTABLE STATEMENT  DQXGE
  1776.       EPMACH = D1MACH(4)
  1777.       UFLOW = D1MACH(1)
  1778. C
  1779. C           TEST ON VALIDITY OF PARAMETERS
  1780. C           ------------------------------
  1781. C
  1782.       IER = 0
  1783.       LAST = 0
  1784.       RESULT = 0.0D+00
  1785.       ABSERR = 0.0D+00
  1786.       ALIST(1) = A
  1787.       BLIST(1) = B
  1788.       RLIST(1) = 0.0D+00
  1789.       ELIST(1) = 0.0D+00
  1790.       IORD(1) = 0
  1791.       IF(EPSABS.LE.0.0D+00.AND.
  1792.      *  EPSREL.LT.DMAX1(0.5D+02*EPMACH,0.5D-28)) IER = 6
  1793.       IF(IER.EQ.6) GO TO 999
  1794. C
  1795. C           FIRST APPROXIMATION TO THE INTEGRAL
  1796. C           -----------------------------------
  1797. C
  1798.       LP(1)=1
  1799.       LN(1)=1
  1800.       VALP(1,1)=F((A+B)*0.5D0)
  1801.       VALN(1,1)=VALP(1,1)
  1802.       CALL DQXLQM(F,A,B,RESULT,ABSERR,DEFABS,RESABS,
  1803.      *           VALP(1,1),VALN(1,1),LP(1),LN(1),KEY)
  1804.       LAST = 1
  1805.       RLIST(1) = RESULT
  1806.       ELIST(1) = ABSERR
  1807.       IORD(1) = 1
  1808. C
  1809. C           TEST ON ACCURACY.
  1810. C
  1811.       ERRBND = DMAX1(EPSABS,EPSREL*DABS(RESULT))
  1812.       IF(ABSERR.LE.0.5D+02*EPMACH*DEFABS.AND.ABSERR.GT.ERRBND) IER = 2
  1813.       IF(LIMIT.EQ.1) IER = 1
  1814.       IF(IER.NE.0.OR.(ABSERR.LE.ERRBND.AND.ABSERR.NE.RESABS)
  1815.      *  .OR.ABSERR.EQ.0.0D+00) GO TO 999
  1816. C
  1817. C           INITIALIZATION
  1818. C           --------------
  1819. C
  1820. C
  1821.       ERRMAX = ABSERR
  1822.       MAXERR = 1
  1823.       AREA = RESULT
  1824.       ERRSUM = ABSERR
  1825.       NRMAX = 1
  1826.       IROFF1 = 0
  1827.       IROFF2 = 0
  1828. C
  1829. C           MAIN DO-LOOP
  1830. C           ------------
  1831. C
  1832.       DO 30 LAST = 2,LIMIT
  1833. C
  1834. C           BISECT THE SUBINTERVAL WITH THE LARGEST ERROR ESTIMATE.
  1835. C
  1836.         A1 = ALIST(MAXERR)
  1837.         B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR))
  1838.         A2 = B1
  1839.         B2 = BLIST(MAXERR)
  1840.         CALL DQXRRD(F,VALN(1,MAXERR),LN(MAXERR),B1,A1,VN1,VP1,LN1,LP1)
  1841.         CALL DQXLQM(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1,VP1,VN1,LP1,LN1,
  1842.      *              KEY)
  1843.         CALL DQXRRD(F,VALP(1,MAXERR),LP(MAXERR),A2,B2,VP2,VN2,LP2,LN2)
  1844.         CALL DQXLQM(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2,VP2,VN2,LP2,LN2,
  1845.      *              KEY)
  1846. C
  1847. C           IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
  1848. C           AND ERROR AND TEST FOR ACCURACY.
  1849. C
  1850.         AREA12 = AREA1+AREA2
  1851.         ERRO12 = ERROR1+ERROR2
  1852.         ERRSUM = ERRSUM+ERRO12-ERRMAX
  1853.         AREA = AREA+AREA12-RLIST(MAXERR)
  1854.         IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2) GO TO 5
  1855.         IF(DABS(RLIST(MAXERR)-AREA12).LE.0.1D-04*DABS(AREA12)
  1856.      *  .AND.ERRO12.GE.0.99D+00*ERRMAX) IROFF1 = IROFF1+1
  1857.         IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF2 = IROFF2+1
  1858.     5   RLIST(MAXERR) = AREA1
  1859.         RLIST(LAST) = AREA2
  1860.         ERRBND = DMAX1(EPSABS,EPSREL*DABS(AREA))
  1861.         IF(ERRSUM.LE.ERRBND) GO TO 8
  1862. C
  1863. C           TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
  1864. C
  1865.         IF(IROFF1.GE.6.OR.IROFF2.GE.20) IER = 2
  1866. C
  1867. C           SET ERROR FLAG IN THE CASE THAT THE NUMBER OF SUBINTERVALS
  1868. C           EQUALS LIMIT.
  1869. C
  1870.         IF(LAST.EQ.LIMIT) IER = 1
  1871. C
  1872. C           SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
  1873. C           AT A POINT OF THE INTEGRATION RANGE.
  1874. C
  1875.         IF(DMAX1(DABS(A1),DABS(B2)).LE.(0.1D+01+0.1D+03*
  1876.      *  EPMACH)*(DABS(A2)+0.1D+04*UFLOW)) IER = 3
  1877. C
  1878. C           APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
  1879. C
  1880.     8   IF(ERROR2.GT.ERROR1) GO TO 10
  1881.         ALIST(LAST) = A2
  1882.         BLIST(MAXERR) = B1
  1883.         BLIST(LAST) = B2
  1884.         ELIST(MAXERR) = ERROR1
  1885.         ELIST(LAST) = ERROR2
  1886.         CALL DQXCPY(VALP(1,MAXERR),VP1,LP1)
  1887.         LP(MAXERR)=LP1
  1888.         CALL DQXCPY(VALN(1,MAXERR),VN1,LN1)
  1889.         LN(MAXERR)=LN1
  1890.         CALL DQXCPY(VALP(1,LAST),VP2,LP2)
  1891.         LP(LAST)=LP2
  1892.         CALL DQXCPY(VALN(1,LAST),VN2,LN2)
  1893.         LN(LAST)=LN2
  1894.         GO TO 20
  1895.    10   ALIST(MAXERR) = A2
  1896.         ALIST(LAST) = A1
  1897.         BLIST(LAST) = B1
  1898.         RLIST(MAXERR) = AREA2
  1899.         RLIST(LAST) = AREA1
  1900.         ELIST(MAXERR) = ERROR2
  1901.         ELIST(LAST) = ERROR1
  1902.         CALL DQXCPY(VALP(1,MAXERR),VP2,LP2)
  1903.         LP(MAXERR)=LP2
  1904.         CALL DQXCPY(VALN(1,MAXERR),VN2,LN2)
  1905.         LN(MAXERR)=LN2
  1906.         CALL DQXCPY(VALP(1,LAST),VP1,LP1)
  1907.         LP(LAST)=LP1
  1908.         CALL DQXCPY(VALN(1,LAST),VN1,LN1)
  1909.         LN(LAST)=LN1
  1910. C
  1911. C           CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING
  1912. C           IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
  1913. C           WITH THE LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
  1914. C
  1915.    20   CALL DQPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX)
  1916. C ***JUMP OUT OF DO-LOOP
  1917.         IF(IER.NE.0.OR.ERRSUM.LE.ERRBND) GO TO 40
  1918.    30 CONTINUE
  1919. C
  1920. C           COMPUTE FINAL RESULT.
  1921. C           ---------------------
  1922. C
  1923.    40 RESULT = 0.0D+00
  1924.       DO 50 K=1,LAST
  1925.         RESULT = RESULT+RLIST(K)
  1926.    50 CONTINUE
  1927.       ABSERR = ERRSUM
  1928.   999 RETURN
  1929.       END
  1930. *DECK DQXLQM
  1931.       SUBROUTINE DQXLQM(F,A,B,RESULT,ABSERR,RESABS,RESASC,VR,VS,LR,LS,
  1932.      *                  KEY)
  1933. C
  1934. C            TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR
  1935. C                           ESTIMATE
  1936. C                       J = INTEGRAL OF ABS(F) OVER (A,B)
  1937. C
  1938. C           PARAMETERS
  1939. C            ON ENTRY
  1940. C              F      - DOUBLE PRECISION
  1941. C                       FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  1942. C                       FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  1943. C                       DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  1944. C
  1945. C              A      - DOUBLE PRECISION
  1946. C                       LOWER LIMIT OF INTEGRATION
  1947. C
  1948. C              B      - DOUBLE PRECISION
  1949. C                       UPPER LIMIT OF INTEGRATION
  1950. C
  1951. C              VR     - DOUBLE PRECISION
  1952. C                       VECTOR OF LENGTH LR CONTAINING THE
  1953. C                       SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS
  1954. C
  1955. C              VS     - DOUBLE PRECISION
  1956. C                       VECTOR OF LENGTH LS CONTAINING THE
  1957. C                       SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS
  1958. C
  1959. C              LR     - INTEGER
  1960. C              LS       NUMBER OF ELEMENTS IN
  1961. C                       VR,VS RESPECTIVELY
  1962. C
  1963. C            KEY    - INTEGER
  1964. C                     KEY FOR CHOICE OF LOCAL INTEGRATION RULE
  1965. C                     RMS FORMULAS ARE USED WITH
  1966. C                      13 - 19               POINTS IF KEY.LT.1,
  1967. C                      13 - 19 - (27)        POINTS IF KEY = 1,
  1968. C                      13 - 19 - (27) - (41) POINTS IF KEY = 2,
  1969. C                           19 -  27  - (41) POINTS IF KEY = 3,
  1970. C                                 27  -  41  POINTS IF KEY.GT.3.
  1971. C
  1972. C                         (RULES) USED IF THE FUNCTION APPEARS
  1973. C                         ENOUGH REGULAR
  1974. C
  1975. C            ON RETURN
  1976. C              RESULT - DOUBLE PRECISION
  1977. C                       APPROXIMATION TO THE INTEGRAL I
  1978. C
  1979. C              ABSERR - DOUBLE PRECISION
  1980. C                       ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
  1981. C                       WHICH SHOULD NOT EXCEED ABS(I-RESULT)
  1982. C
  1983. C              RESABS - DOUBLE PRECISION
  1984. C                       APPROXIMATION TO THE INTEGRAL J
  1985. C
  1986. C              RESASC - DOUBLE PRECISION
  1987. C                       APPROXIMATION TO THE INTEGRAL OF ABS(F-I/(B-A))
  1988. C                       OVER (A,B)
  1989. C
  1990. C              VR     - DOUBLE PRECISION
  1991. C                       VECTOR OF LENGTH LR CONTAINING THE
  1992. C                       SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS
  1993. C
  1994. C              VS     - DOUBLE PRECISION
  1995. C                       VECTOR OF LENGTH LS CONTAINING THE
  1996. C                       SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS
  1997. C
  1998. C              LR     - INTEGER
  1999. C              LS       NUMBER OF ELEMENTS IN
  2000. C                       VR,VS RESPECTIVELY
  2001. C
  2002. C***ROUTINES CALLED  D1MACH,DQXRUL
  2003. C
  2004.       DOUBLE PRECISION F,A,B,RESULT,ABSERR,RESABS,RESASC,
  2005.      *          D1MACH,EPMACH,RESG,RESK,UFLOW,ERROLD,VR(21),VS(21)
  2006.       INTEGER K,K0,K1,K2,KEY,KEY1,LR,LS
  2007.       EXTERNAL F
  2008. C
  2009. C            MACHINE DEPENDENT CONSTANTS
  2010. C            ---------------------------
  2011. C
  2012. C          EPMACH IS THE LARGEST RELATIVE SPACING.
  2013. C          UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
  2014. C          ERROLD IS THE LARGEST MAGNITUDE.
  2015. C
  2016. C***FIRST EXECUTABLE STATEMENT DQXLQM
  2017.       EPMACH = D1MACH(4)
  2018.       UFLOW  = D1MACH(1)
  2019.       ERROLD = D1MACH(2)
  2020. C
  2021.       KEY1 = MAX(KEY ,  0)
  2022.       KEY1 = MIN(KEY1,  4)
  2023.       K0   = MAX(KEY1-2,0)
  2024.       K1   = K0+1
  2025.       K2   = MIN(KEY1+1,3)
  2026. C
  2027.       CALL DQXRUL(F,A,B,RESG,RESABS,RESASC,K0,K1,VR,VS,LR,LS)
  2028.       DO 99 K=K1,K2
  2029.         CALL DQXRUL(F,A,B,RESK,RESABS,RESASC,K,K1,VR,VS,LR,LS)
  2030.         RESULT=RESK
  2031.         ABSERR = DABS((RESK-RESG))
  2032.         IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00)
  2033.      *  ABSERR = RESASC*DMIN1(1.D0,(200D0*ABSERR/RESASC)**1.5D+00)
  2034.         IF(RESABS.GT.UFLOW/(10D0*EPMACH)) ABSERR = DMAX1
  2035.      *  ((EPMACH*10D0)*RESABS,ABSERR)
  2036.         RESG=RESK
  2037.         IF(ABSERR.GT.ERROLD*0.16)GOTO 3000
  2038.         IF(ABSERR.LT.1000*EPMACH*RESABS)GOTO 3000
  2039.         ERROLD=ABSERR
  2040. 99    CONTINUE
  2041. 3000  CONTINUE
  2042.       RETURN
  2043.       END
  2044. *DECK DQXRUL
  2045.       SUBROUTINE DQXRUL(F,XL,XU,Y,YA,YM,KE,K1,FV1,FV2,L1,L2)
  2046. C
  2047. C            TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR
  2048. C                           ESTIMATE
  2049. C            AND CONDITIONALLY COMPUTE
  2050. C                       J = INTEGRAL OF ABS(F) OVER (A,B)
  2051. C                       BY USING AN  RMS RULE
  2052. C           PARAMETERS
  2053. C            ON ENTRY
  2054. C              F      - DOUBLE PRECISION
  2055. C                       FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  2056. C                       FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  2057. C                       DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  2058. C
  2059. C              XL     - DOUBLE PRECISION
  2060. C                       LOWER LIMIT OF INTEGRATION
  2061. C
  2062. C              XU     - DOUBLE PRECISION
  2063. C                       UPPER LIMIT OF INTEGRATION
  2064. C
  2065. C
  2066. C              KE     - INTEGER
  2067. C                     KEY FOR CHOICE OF LOCAL INTEGRATION RULE
  2068. C                     AN RMS RULE IS USED WITH
  2069. C                         13      POINTS IF KE  = 2,
  2070. C                         19      POINTS IF KE  = 3,
  2071. C                         27      POINTS IF KE  = 4,
  2072. C                         42      POINTS IF KE  = 5
  2073. C
  2074. C              K1     INTEGER
  2075. C                     VALUE OF KEY FOR WHICH THE ADDITIONAL ESTIMATES
  2076. C                     YA, YM ARE TO BE COMPUTED
  2077. C
  2078. C              FV1    - DOUBLE PRECISION
  2079. C                       VECTOR CONTAINING L1
  2080. C                       SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS
  2081. C
  2082. C              FV2    - DOUBLE PRECISION
  2083. C                       VECTOR CONTAINING L2
  2084. C                       SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS
  2085. C
  2086. C              L1     - INTEGER
  2087. C              L2       NUMBER OF ELEMENTS IN FV1,FV2  RESPECTIVELY
  2088. C
  2089. C            ON RETURN
  2090. C              Y      - DOUBLE PRECISION
  2091. C                       APPROXIMATION TO THE INTEGRAL I
  2092. C                       RESULT IS COMPUTED BY APPLYING THE
  2093. C                       REQUESTED RMS RULE
  2094. C
  2095. C              YA     - DOUBLE PRECISION
  2096. C                       IF KEY = K1  APPROXIMATION TO THE INTEGRAL J
  2097. C                       ELSE UNCHANGED
  2098. C
  2099. C              YM     - DOUBLE PRECISION
  2100. C                       IF KEY = K1  APPROXIMATION TO THE INTEGRAL OF
  2101. C                                      ABS(F-I/(XU-XL)   OVER (XL,XU)
  2102. C                       ELSE UNCHANGED
  2103. C
  2104. C              FV1    - DOUBLE PRECISION
  2105. C                       VECTOR L1 CONTAINING L1
  2106. C                       SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS
  2107. C
  2108. C              FV2    - DOUBLE PRECISION
  2109. C                       VECTOR CONTAINING L2
  2110. C                       SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS
  2111. C
  2112. C              L1     - INTEGER
  2113. C              L2       NUMBER OF ELEMENTS IN FV1,FV2  RESPECTIVELY
  2114. C
  2115. C***ROUTINES CALLED  F
  2116. C
  2117.       DOUBLE PRECISION F,XL,XU,LDL,Y,YA,YM,Y2,XX(41),WW(52),
  2118.      *                 FV1(21),FV2(21),AA,BB,C
  2119.       EXTERNAL F
  2120.       INTEGER ISTART(4),LEN(4),KE,K1,L1,L2
  2121.       SAVE ISTART,LEN,XX,WW
  2122.       DATA ISTART/0, 7, 17, 31/,LEN/7, 10, 14, 21/
  2123.       DATA XX(  1)/.0                       /
  2124.       DATA XX(  2)/.25000000000000000000D+00/
  2125.       DATA XX(  3)/.50000000000000000000D+00/
  2126.       DATA XX(  4)/.75000000000000000000D+00/
  2127.       DATA XX(  5)/.87500000000000000000D+00/
  2128.       DATA XX(  6)/.93750000000000000000D+00/
  2129.       DATA XX(  7)/.10000000000000000000D+01/
  2130.       DATA XX(  8)/.37500000000000000000D+00/
  2131.       DATA XX(  9)/.62500000000000000000D+00/
  2132.       DATA XX( 10)/.96875000000000000000D+00/
  2133.       DATA XX( 11)/.12500000000000000000D+00/
  2134.       DATA XX( 12)/.68750000000000000000D+00/
  2135.       DATA XX( 13)/.81250000000000000000D+00/
  2136.       DATA XX( 14)/.98437500000000000000D+00/
  2137.       DATA XX( 15)/.18750000000000000000D+00/
  2138.       DATA XX( 16)/.31250000000000000000D+00/
  2139.       DATA XX( 17)/.43750000000000000000D+00/
  2140.       DATA XX( 18)/.56250000000000000000D+00/
  2141.       DATA XX( 19)/.84375000000000000000D+00/
  2142.       DATA XX( 20)/.90625000000000000000D+00/
  2143.       DATA XX( 21)/.99218750000000000000D+00/
  2144. C   NUMBER OF NODES 13
  2145.       DATA WW(1)/1.303262173284849021810473057638590518409112513421D-1/
  2146.       DATA WW(2)/2.390632866847646220320329836544615917290026806242D-1/
  2147.       DATA WW(3)/2.630626354774670227333506083741355715758124943143D-1/
  2148.       DATA WW(4)/2.186819313830574175167853094864355208948886875898D-1/
  2149.       DATA WW(5)/2.757897646642836865859601197607471574336674206700D-2/
  2150.       DATA WW(6)/1.055750100538458443365034879086669791305550493830D-1/
  2151.       DATA WW(7)/1.571194260595182254168429283636656908546309467968D-2/
  2152. C   NUMBER OF NODES 19
  2153.       DATA WW(8)/1.298751627936015783241173611320651866834051160074D-1/
  2154.       DATA WW(9)/2.249996826462523640447834514709508786970828213187D-1/
  2155.       DATA WW(15)/5.542699233295875168406783695143646338274805359780D-2/
  2156.       DATA WW(10)/1.680415725925575286319046726692683040162290325505D-1/
  2157.       DATA WW(16)/9.986735247403367525720377847755415293097913496236D-2/
  2158.       DATA WW(11)/1.415567675701225879892811622832845252125600939627D-1/
  2159.       DATA WW(12)/1.006482260551160175038684459742336605269707889822D-1/
  2160.       DATA WW(13)/2.510604860724282479058338820428989444699235030871D-2/
  2161.       DATA WW(17)/4.507523056810492466415880450799432587809828791196D-2/
  2162.       DATA WW(14)/9.402964360009747110031098328922608224934320397592D-3/
  2163. C   NUMBER OF NODES 27
  2164.       DATA WW(18)/6.300942249647773931746170540321811473310938661469D-2/
  2165.       DATA WW(28)/1.239572396231834242194189674243818619042280816640D-1/
  2166.       DATA WW(19)/1.261383225537664703012999637242003647020326905948D-1/
  2167.       DATA WW(25)/1.235837891364555000245004813294817451524633100256D-1/
  2168.       DATA WW(20)/1.273864433581028272878709981850307363453523117880D-1/
  2169.       DATA WW(26)/1.148933497158144016800199601785309838604146040215D-1/
  2170.       DATA WW(29)/2.501306413750310579525950767549691151739047969345D-2/
  2171.       DATA WW(21)/8.576500414311820514214087864326799153427368592787D-2/
  2172.       DATA WW(30)/4.915957918146130094258849161350510503556792927578D-2/
  2173.       DATA WW(22)/7.102884842310253397447305465997026228407227220665D-2/
  2174.       DATA WW(23)/5.026383572857942403759829860675892897279675661654D-2/
  2175.       DATA WW(27)/1.252575774226122633391477702593585307254527198070D-2/
  2176.       DATA WW(31)/2.259167374956474713302030584548274729936249753832D-2/
  2177.       DATA WW(24)/4.683670010609093810432609684738393586390722052124D-3/
  2178. C   NUMBER OF NODES 41
  2179.       DATA WW(32)/6.362762978782724559269342300509058175967124446839D-2/
  2180.       DATA WW(42)/1.187141856692283347609436153545356484256869129472D-1/
  2181.       DATA WW(46)/1.533126874056586959338368742803997744815413565014D-2/
  2182.       DATA WW(33)/9.950065827346794643193261975720606296171462239514D-2/
  2183.       DATA WW(47)/3.527159369750123100455704702965541866345781113903D-2/
  2184.       DATA WW(39)/8.140326425945938045967829319725797511040878579808D-2/
  2185.       DATA WW(48)/5.000556431653955124212795201196389006184693561679D-2/
  2186.       DATA WW(34)/7.048220002718565366098742295389607994441704889441D-2/
  2187.       DATA WW(49)/5.744164831179720106340717579281831675999717767532D-2/
  2188.       DATA WW(40)/6.583213447600552906273539578430361199084485578379D-2/
  2189.       DATA WW(43)/5.999947605385971985589674757013565610751028128731D-2/
  2190.       DATA WW(35)/6.512297339398335645872697307762912795346716454337D-2/
  2191.       DATA WW(44)/5.500937980198041736910257988346101839062581489820D-2/
  2192.       DATA WW(50)/1.598823797283813438301248206397233634639162043386D-2/
  2193.       DATA WW(36)/3.998229150313659724790527138690215186863915308702D-2/
  2194.       DATA WW(51)/2.635660410220884993472478832884065450876913559421D-2/
  2195.       DATA WW(37)/3.456512257080287509832054272964315588028252136044D-2/
  2196.       DATA WW(41)/2.592913726450792546064232192976262988065252032902D-2/
  2197.       DATA WW(45)/5.264422421764655969760271538981443718440340270116D-3/
  2198.       DATA WW(52)/1.196003937945541091670106760660561117114584656319D-2/
  2199.       DATA WW(38)/2.212167975884114432760321569298651047876071264944D-3/
  2200. C
  2201. C***FIRST EXECUTABLE STATEMENT DQXRUL
  2202.       K=KE+1
  2203.       IS=ISTART(K)
  2204.       KS=LEN(K)
  2205.       LDL= XU-XL
  2206.       BB= LDL*0.5D0
  2207.       AA= XL+BB
  2208.       Y =0.D0
  2209.       DO 100 I=1,KS
  2210.          IF(I.GT.L1.OR.I.GT.L2)   C=BB*XX(I)
  2211.          IF(I.GT.L1)              FV1(I)=F(AA+C)
  2212.          IF(I.GT.L2)              FV2(I)=F(AA-C)
  2213. 100      Y=Y+(FV1(I)+FV2(I))*WW(IS+I)
  2214.       Y2=Y
  2215.       Y=Y*BB
  2216.       IF(L1.LT.KS)L1=KS
  2217.       IF(L2.LT.KS)L2=KS
  2218.       IF(KE.NE.K1)RETURN
  2219.       YA=0.D0
  2220.       DO 25 I=1,KS
  2221.   25    YA=YA+(DABS(FV1(I))+DABS(FV2(I)))*WW(IS+I)
  2222.       Y2=Y2*0.5D0
  2223.       YM=0.D0
  2224.       YA=YA*DABS(BB)
  2225.       DO 27 I=1,KS
  2226.    27 YM=YM+(DABS(FV1(I)-Y2)+DABS(FV2(I)-Y2))*WW(IS+I)
  2227.       YM=YM*DABS(BB)
  2228.       RETURN
  2229.       END
  2230. *DECK DQXRRD
  2231.       SUBROUTINE DQXRRD(F,Z,LZ,XL,XU,R,S,LR,LS)
  2232. C
  2233. C            TO REORDER THE COMPUTED FUNCTIONAL VALUES BEFORE
  2234. C            THE BISECTION OF AN INTERVAL
  2235. C           PARAMETERS
  2236. C            ON ENTRY
  2237. C              F      - DOUBLE PRECISION
  2238. C                       FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  2239. C                       FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  2240. C                       DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  2241. C
  2242. C              XL     - DOUBLE PRECISION
  2243. C                       LOWER LIMIT OF INTERVAL
  2244. C
  2245. C              XU     - DOUBLE PRECISION
  2246. C                       UPPER LIMIT OF INTERVAL
  2247. C
  2248. C              Z      - DOUBLE PRECISION
  2249. C                       VECTOR CONTAINING LZ
  2250. C                       SAVED  FUNCTIONAL VALUES
  2251. C
  2252. C              LZ     - INTEGER
  2253. C                       NUMBER OF ELEMENTS IN LZ
  2254. C
  2255. C
  2256. C            ON RETURN
  2257. C              R      - DOUBLE PRECISION
  2258. C              S        VECTORS CONTAINING LR, LS
  2259. C                       SAVED  FUNCTIONAL VALUES FOR THE NEW INTERVALS
  2260. C
  2261. C              LR     - INTEGER
  2262. C              LS       NUMBER OF ELEMENTES IN R,S RESPECTIVELY
  2263. C
  2264. C***ROUTINES CALLED  F
  2265. C
  2266.       DOUBLE PRECISION F,R,S,Z,XU,XL,DLEN,CENTR
  2267.       INTEGER LR,LS,LZ
  2268.       DIMENSION R(21),S(21),Z(21)
  2269. C***FIRST EXECUTABLE STATEMENT DQXRRD
  2270.       DLEN= (XU-XL)*0.5D0
  2271.       CENTR= XL+DLEN
  2272.       R(1)=  Z(3)
  2273.       R(2)=  Z(9)
  2274.       R(3)=  Z(4)
  2275.       R(4)=  Z(5)
  2276.       R(5)=  Z(6)
  2277.       R(6)=  Z(10)
  2278.       R(7)=  Z(7)
  2279.       S(1)=  Z(3)
  2280.       S(2)=  Z(8)
  2281.       S(3)=  Z(2)
  2282.       S(7)=  Z(1)
  2283.       IF(LZ.GT.11)GOTO 2
  2284.       R(8)=  F(CENTR+DLEN*0.37500000D0)
  2285.       R(9)=  F(CENTR+DLEN*0.62500000D0)
  2286.       R(10)=  F(CENTR+DLEN*0.96875000D0)
  2287.       LR=  10
  2288.       IF(LZ.NE.11)S(4)=  F(CENTR-DLEN*0.75000000D0)
  2289.       IF(LZ.EQ.11)S(4)=  Z(11)
  2290.       S(5)=  F(CENTR-DLEN*0.87500000D0)
  2291.       S(6)=  F(CENTR-DLEN*0.93750000D0)
  2292.       S(8)=  F(CENTR-DLEN*0.37500000D0)
  2293.       S(9)=  F(CENTR-DLEN*0.62500000D0)
  2294.       S(10)=  F(CENTR-DLEN*0.96875000D0)
  2295.       LS=  10
  2296.       GOTO 3000
  2297. 2     R(8)= Z(12)
  2298.       R(9)= Z(13)
  2299.       R(10)= Z(14)
  2300.       LR=  10
  2301.       S(4)= Z(11)
  2302.       S(5)= F(CENTR-DLEN*0.87500000D0)
  2303.       S(6)= F(CENTR-DLEN*0.93750000D0)
  2304.       IF(LZ.GT.14)GOTO3
  2305.       S(8)= F(CENTR-DLEN*0.37500000D0)
  2306.       S(9)= F(CENTR-DLEN*0.62500000D0)
  2307.       S(10)= F(CENTR-DLEN*0.96875000D0)
  2308.       LS=  10
  2309.       GOTO 3000
  2310. 3     R(11)= Z(18)
  2311.       R(12)= Z(19)
  2312.       R(13)= Z(20)
  2313.       R(14)= Z(21)
  2314.       LR=  14
  2315.       S(8)= Z(16)
  2316.       S(9)= Z(15)
  2317.       S(10)= F(CENTR-DLEN*0.96875000D0)
  2318.       S(11)= Z(17)
  2319.       LS=  11
  2320. 3000  RETURN
  2321.       END
  2322. *DECK DQXCPY
  2323.       SUBROUTINE DQXCPY(A,B,L)
  2324. C
  2325. C  TO COPY THE DOUBLE PRECISION VECTOR B OF LENGTH L   I N T O
  2326. C          THE DOUBLE PRECISION VECTOR A OF LENGTH L
  2327. C
  2328. C***REMARK  THIS ROUTINE CAN BE IMPROVED, BY CODING IT IN THE
  2329. C           ASSEMBLER LANGUAGE OF THE USED MACHINE
  2330. C***ROUTINES CALLED  (NONE)
  2331. C
  2332.       INTEGER L
  2333.       DOUBLE PRECISION A(L),B(L)
  2334. C***FIRST EXECUTABLE STATEMENT DQXCPY
  2335.       DO 1 I=1,L
  2336. 1     A(I)=B(I)
  2337.       RETURN
  2338.       END
  2339. *DECK DTEST  *** TEST PROGRAM ***
  2340. C  THIS DRIVER CALLS VARIOUS INTEGRATORS WITH  VARIOUS RELATIVE
  2341. C  TOLERACES TO INTEGRATE THE FUNCTION
  2342. C     F(X) = DSQRT(X/(1-X)) * DLOG(X)   OVER (0, 1)
  2343. C
  2344. C THE EXACT INTEGRAL IS                           -0.606789763508705...
  2345. C
  2346.       DOUBLE PRECISION F,A,B,ATOL,RTOL,TEE,Y,WORK(4600)
  2347.       INTEGER NCALLS,IER,NEVAL,LAST,IWORK(300),KKK,KEY,LIMIT,LENW,LENIW
  2348.       EXTERNAL F
  2349.       COMMON /CTEST/NCALLS
  2350.       DATA LIMIT/100/,LENIW/300/,LENW/4600/
  2351.       WRITE(6,600)
  2352. 600   FORMAT(32H1 COMPUTATION OF THE INTEGRAL OF,
  2353.      *      /38H       F(X) = DSQRT(X/(1-X)) * DLOG(X),
  2354.      *      /13H  OVER (0, 1),
  2355.      *      /23H  THE EXACT INTEGRAL IS,30X,
  2356.      *       22H -0.606789763508705...)
  2357.       A=0.0D0
  2358.       B=1.0D0
  2359.       WRITE(6,61)
  2360.  61   FORMAT(//22H TEST OF DQXGS               //
  2361.      *  53H RELATIVE TOLERANCE FUNCTION CALLS ERROR ESTIMATE IER,
  2362.      *  13H       RESULT )
  2363.       DO 1 I=3,13,2
  2364.          ATOL=1D-20
  2365.          RTOL=10D0**(-I)
  2366.          NCALLS=0
  2367.          CALL DQXGS (F,A,B,ATOL,RTOL,Y,TEE,IER,LIMIT,
  2368.      *               LENIW,LENW,LAST,IWORK,WORK)
  2369.          WRITE(6,60)I,NCALLS,TEE,IER,Y
  2370. 60       FORMAT(4X,6H10**(-,I2,1H),I17,9X,D7.2,4X,I2,F20.15)
  2371.  1    CONTINUE
  2372.       DO 2 KKK = 1,5
  2373.        KEY = KKK - 1
  2374.        WRITE(6,62) KEY
  2375.  62    FORMAT(//23H TEST OF DQXG , KEY =    ,I4   //
  2376.      *  53H RELATIVE TOLERANCE FUNCTION CALLS ERROR ESTIMATE IER,
  2377.      *  13H       RESULT )
  2378.        DO 2 I=3,13,2
  2379.          ATOL=1D-20
  2380.          RTOL=10D0**(-I)
  2381.          NCALLS=0
  2382.          CALL DQXG  (F,A,B,ATOL,RTOL,KEY,Y,TEE,IER,LIMIT,
  2383.      *               LENIW,LENW,LAST,IWORK,WORK)
  2384.          WRITE(6,60)I,NCALLS,TEE,IER,Y
  2385.  2    CONTINUE
  2386.       STOP
  2387.       END
  2388. *DECK F
  2389.       DOUBLE PRECISION FUNCTION F(X)
  2390.       DOUBLE PRECISION X
  2391. C        INTEGRAND FUNCTION, IN THE INTERVAL (0,1) THE EXACT INTEGRAL
  2392. C   TO 15 DIGITS IS
  2393. C          -0.60678 97635 08705D0
  2394. C
  2395.       COMMON /CTEST/NCALLS
  2396.       NCALLS=NCALLS+1
  2397.       IF (X.EQ.0D0) GO TO  10
  2398.       IF (X.EQ.1D0) GO TO  10
  2399.       F=DSQRT(X/(1-X))*DLOG(X)
  2400.       RETURN
  2401. 10    F=0.0
  2402.       RETURN
  2403.       END
  2404. CC+------------------------------------------------------------------+CC
  2405. CC+                                                                  +CC
  2406. CC+                LAWRENCE LIVERMORE NATIONAL LABORATORY            +CC
  2407. CC+                                                                  +CC
  2408. CC+        MATHEMATICS AND STATISTICS DIVISION SOFTWARE LIBRARY      +CC
  2409. CC+                                                                  +CC
  2410. CC+------------------------------------------------------------------+CC
  2411. CC+                                                                  +CC
  2412. CC+  THE SLATEC MATHEMATICAL PROGRAM LIBRARY IS ISSUED BY            +CC
  2413. CC+  THE FOLLOWING:                                                  +CC
  2414. CC+                                                                  +CC
  2415. CC+       AIR FORCE WEAPONS LABORATORY, ALBUQUERQUE                  +CC
  2416. CC+       LAWRENCE LIVERMORE NATIONAL LABORATORY, LIVERMORE          +CC
  2417. CC+       LOS ALAMOS NATIONAL SCIENTIFIC LABORATORY, LOS ALAMOS      +CC
  2418. CC+       NATIONAL BUREAU OF STANDARDS, WASHINGTON                   +CC
  2419. CC+       OAK RIDGE NATIONAL LABORATORY, OAK RIDGE                   +CC
  2420. CC+       SANDIA NATIONAL LABORATORIES, ALBUQUERQUE                  +CC
  2421. CC+       SANDIA NATIONAL LABORATORIES, LIVERMORE                    +CC
  2422. CC+                                                                  +CC
  2423. CC+  PLEASE REFER QUESTIONS REGARDING THIS SOFTWARE TO THE           +CC
  2424. CC+  MSD CONSULTING OFFICE, EXTENSION 3-2976.                        +CC
  2425. CC+                                                                  +CC
  2426. CC+                                                                  +CC
  2427. CC+              * * * * * N O T I C E * * * * *                     +CC
  2428. CC+                                                                  +CC
  2429. CC+  THIS MATERIAL WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED      +CC
  2430. CC+  BY THE UNITED STATES GOVERNMENT.  NEITHER THE UNITED            +CC
  2431. CC+  STATES GOVERNMENT, NOR THE DEPARTMENT OF DEFENSE, NOR           +CC
  2432. CC+  ANY OF THEIR EMPLOYEES, MAKES ANY WARRANTY, EXPRESSED           +CC
  2433. CC+  OR IMPLIED, OR ASSUMES ANY LEGAL LIABILITY OR                   +CC
  2434. CC+  RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR               +CC
  2435. CC+  USEFULNESS OF ANY INFORMATION, APPARATUS, PRODUCT,              +CC
  2436. CC+  OR PROCESS DISCLOSED, OR REPRESENTS ITS USE WOULD NOT           +CC
  2437. CC+  INFRINGE UPON PRIVATELY OWNED RIGHTS.                           +CC
  2438. CC+                                                                  +CC
  2439. CC+------------------------------------------------------------------+CC
  2440. *DECK QPSRT
  2441.       SUBROUTINE QPSRT(LIMIT,LAST,MAXERR,ERMAX,ELIST,IORD,NRMAX)
  2442. C***BEGIN PROLOGUE  QPSRT
  2443. C***REFER TO  QAGE,QAGIE,QAGPE,QAGSE,QAWCE,QAWOE,QAWSE
  2444. C***ROUTINES CALLED  (NONE)
  2445. C***KEYWORDS  SEQUENTIAL SORTING
  2446. C***DESCRIPTION
  2447. C
  2448. C 1.        QPSRT
  2449. C           ORDERING ROUTINE
  2450. C              STANDARD FORTRAN SUBROUTINE
  2451. C              REAL VERSION
  2452. C
  2453. C 2.        PURPOSE
  2454. C              THIS ROUTINE MAINTAINS THE DESCENDING ORDERING
  2455. C              IN THE LIST OF THE LOCAL ERROR ESTIMATES RESULTING FROM
  2456. C              THE INTERVAL SUBDIVISION PROCESS. AT EACH CALL TWO ERROR
  2457. C              ESTIMATES ARE INSERTED USING THE SEQUENTIAL SEARCH
  2458. C              METHOD, TOP-DOWN FOR THE LARGEST ERROR ESTIMATE
  2459. C              AND BOTTOM-UP FOR THE SMALLEST ERROR ESTIMATE.
  2460. C
  2461. C 3.        CALLING SEQUENCE
  2462. C              CALL QPSRT(LIMIT,LAST,MAXERR,ERMAX,ELIST,IORD,NRMAX)
  2463. C
  2464. C           PARAMETERS (MEANING AT OUTPUT)
  2465. C              LIMIT  - INTEGER
  2466. C                       MAXIMUM NUMBER OF ERROR ESTIMATES THE LIST
  2467. C                       CAN CONTAIN
  2468. C
  2469. C              LAST   - INTEGER
  2470. C                       NUMBER OF ERROR ESTIMATES CURRENTLY
  2471. C                       IN THE LIST
  2472. C
  2473. C              MAXERR - INTEGER
  2474. C                       MAXERR POINTS TO THE NRMAX-TH LARGEST ERROR
  2475. C                       ESTIMATE CURRENTLY IN THE LIST
  2476. C
  2477. C              ERMAX  - REAL
  2478. C                       NRMAX-TH LARGEST ERROR ESTIMATE
  2479. C                       ERMAX = ELIST(MAXERR)
  2480. C
  2481. C              ELIST  - REAL
  2482. C                       VECTOR OF DIMENSION LAST CONTAINING
  2483. C                       THE ERROR ESTIMATES
  2484. C
  2485. C              IORD   - INTEGER
  2486. C                       VECTOR OF DIMENSION LAST, THE FIRST K
  2487. C                       ELEMENTS OF WHICH CONTAIN POINTERS
  2488. C                       TO THE ERROR ESTIMATES, SUCH THAT
  2489. C                       ELIST(IORD(1)),... , ELIST(IORD(K))
  2490. C                       FORM A DECREASING SEQUENCE, WITH
  2491. C                       K = LAST IF LAST.LE.(LIMIT/2+2), AND
  2492. C                       K = LIMIT+1-LAST OTHERWISE
  2493. C
  2494. C              NRMAX  - INTEGER
  2495. C                       MAXERR = IORD(NRMAX)
  2496. C
  2497. C 4.        NO SUBROUTINES OR FUNCTIONS NEEDED
  2498. C***END PROLOGUE  QPSRT
  2499. C
  2500.       REAL ELIST,ERMAX,ERRMAX,ERRMIN
  2501.       INTEGER I,IBEG,IDO,IORD,ISUCC,J,JBND,JUPBN,K,LAST,LIMIT,MAXERR,
  2502.      1  NRMAX
  2503.       DIMENSION ELIST(LAST),IORD(LAST)
  2504. C
  2505. C           CHECK WHETHER THE LIST CONTAINS MORE THAN
  2506. C           TWO ERROR ESTIMATES.
  2507. C
  2508. C***FIRST EXECUTABLE STATEMENT  QPSRT
  2509.       IF(LAST.GT.2) GO TO 10
  2510.       IORD(1) = 1
  2511.       IORD(2) = 2
  2512.       GO TO 90
  2513. C
  2514. C           THIS PART OF THE ROUTINE IS ONLY EXECUTED
  2515. C           IF, DUE TO A DIFFICULT INTEGRAND, SUBDIVISION
  2516. C           INCREASED THE ERROR ESTIMATE. IN THE NORMAL CASE
  2517. C           THE INSERT PROCEDURE SHOULD START AFTER THE
  2518. C           NRMAX-TH LARGEST ERROR ESTIMATE.
  2519. C
  2520.    10 ERRMAX = ELIST(MAXERR)
  2521.       IF(NRMAX.EQ.1) GO TO 30
  2522.       IDO = NRMAX-1
  2523.       DO 20 I = 1,IDO
  2524.         ISUCC = IORD(NRMAX-1)
  2525. C ***JUMP OUT OF DO-LOOP
  2526.         IF(ERRMAX.LE.ELIST(ISUCC)) GO TO 30
  2527.         IORD(NRMAX) = ISUCC
  2528.         NRMAX = NRMAX-1
  2529.    20    CONTINUE
  2530. C
  2531. C           COMPUTE THE NUMBER OF ELEMENTS IN THE LIST TO
  2532. C           BE MAINTAINED IN DESCENDING ORDER. THIS NUMBER
  2533. C           DEPENDS ON THE NUMBER OF SUBDIVISIONS STILL
  2534. C           ALLOWED.
  2535. C
  2536.    30 JUPBN = LAST
  2537.       IF(LAST.GT.(LIMIT/2+2)) JUPBN = LIMIT+3-LAST
  2538.       ERRMIN = ELIST(LAST)
  2539. C
  2540. C           INSERT ERRMAX BY TRAVERSING THE LIST TOP-DOWN,
  2541. C           STARTING COMPARISON FROM THE ELEMENT ELIST(IORD(NRMAX+1)).
  2542. C
  2543.       JBND = JUPBN-1
  2544.       IBEG = NRMAX+1
  2545.       IF(IBEG.GT.JBND) GO TO 50
  2546.       DO 40 I=IBEG,JBND
  2547.         ISUCC = IORD(I)
  2548. C ***JUMP OUT OF DO-LOOP
  2549.         IF(ERRMAX.GE.ELIST(ISUCC)) GO TO 60
  2550.         IORD(I-1) = ISUCC
  2551.    40 CONTINUE
  2552.    50 IORD(JBND) = MAXERR
  2553.       IORD(JUPBN) = LAST
  2554.       GO TO 90
  2555. C
  2556. C           INSERT ERRMIN BY TRAVERSING THE LIST BOTTOM-UP.
  2557. C
  2558.    60 IORD(I-1) = MAXERR
  2559.       K = JBND
  2560.       DO 70 J=I,JBND
  2561.         ISUCC = IORD(K)
  2562. C ***JUMP OUT OF DO-LOOP
  2563.         IF(ERRMIN.LT.ELIST(ISUCC)) GO TO 80
  2564.         IORD(K+1) = ISUCC
  2565.         K = K-1
  2566.    70 CONTINUE
  2567.       IORD(I) = LAST
  2568.       GO TO 90
  2569.    80 IORD(K+1) = LAST
  2570. C
  2571. C           SET MAXERR AND ERMAX.
  2572. C
  2573.    90 MAXERR = IORD(NRMAX)
  2574.       ERMAX = ELIST(MAXERR)
  2575.       RETURN
  2576.       END
  2577. *DECK QELG
  2578.       SUBROUTINE QELG(N,EPSTAB,RESULT,ABSERR,RES3LA,NRES)
  2579. C***BEGIN PROLOGUE  QELG
  2580. C***REFER TO  QAGIE,QAGOE,QAGPE,QAGSE
  2581. C***ROUTINES CALLED  R1MACH
  2582. C***REVISION DATE  830518   (YYMMDD)
  2583. C***KEYWORDS  CONVERGENCE ACCELERATION,EPSILON ALGORITHM,EXTRAPOLATION
  2584. C***AUTHOR  PIESSENS, ROBERT, APPLIED MATH. AND PROGR. DIV. -
  2585. C             K. U. LEUVEN
  2586. C           DE DONCKER, ELISE, APPLIED MATH. AND PROGR. DIV. -
  2587. C             K. U. LEUVEN
  2588. C***PURPOSE  THE ROUTINE DETERMINES THE LIMIT OF A GIVEN SEQUENCE OF
  2589. C            APPROXIMATIONS, BY MEANS OF THE EPSILON ALGORITHM OF
  2590. C            P. WYNN. AN ESTIMATE OF THE ABSOLUTE ERROR IS ALSO GIVEN.
  2591. C            THE CONDENSED EPSILON TABLE IS COMPUTED. ONLY THOSE
  2592. C            ELEMENTS NEEDED FOR THE COMPUTATION OF THE NEXT DIAGONAL
  2593. C            ARE PRESERVED.
  2594. C***DESCRIPTION
  2595. C
  2596. C           EPSILON ALGORITHM
  2597. C           STANDARD FORTRAN SUBROUTINE
  2598. C           REAL VERSION
  2599. C
  2600. C           PARAMETERS
  2601. C              N      - INTEGER
  2602. C                       EPSTAB(N) CONTAINS THE NEW ELEMENT IN THE
  2603. C                       FIRST COLUMN OF THE EPSILON TABLE.
  2604. C
  2605. C              EPSTAB - REAL
  2606. C                       VECTOR OF DIMENSION 52 CONTAINING THE ELEMENTS
  2607. C                       OF THE TWO LOWER DIAGONALS OF THE TRIANGULAR
  2608. C                       EPSILON TABLE. THE ELEMENTS ARE NUMBERED
  2609. C                       STARTING AT THE RIGHT-HAND CORNER OF THE
  2610. C                       TRIANGLE.
  2611. C
  2612. C              RESULT - REAL
  2613. C                       RESULTING APPROXIMATION TO THE INTEGRAL
  2614. C
  2615. C              ABSERR - REAL
  2616. C                       ESTIMATE OF THE ABSOLUTE ERROR COMPUTED FROM
  2617. C                       RESULT AND THE 3 PREVIOUS RESULTS
  2618. C
  2619. C              RES3LA - REAL
  2620. C                       VECTOR OF DIMENSION 3 CONTAINING THE LAST 3
  2621. C                       RESULTS
  2622. C
  2623. C              NRES   - INTEGER
  2624. C                       NUMBER OF CALLS TO THE ROUTINE
  2625. C                       (SHOULD BE ZERO AT FIRST CALL)
  2626. C***END PROLOGUE  QELG
  2627. C
  2628.       REAL ABSERR,DELTA1,DELTA2,DELTA3,R1MACH,
  2629.      1  EPMACH,EPSINF,EPSTAB,ERROR,ERR1,ERR2,ERR3,E0,E1,E1ABS,E2,E3,
  2630.      2  OFLOW,RES,RESULT,RES3LA,SS,TOL1,TOL2,TOL3
  2631.       INTEGER I,IB,IB2,IE,INDX,K1,K2,K3,LIMEXP,N,NEWELM,NRES,NUM
  2632.       DIMENSION EPSTAB(52),RES3LA(3)
  2633. C
  2634. C           LIST OF MAJOR VARIABLES
  2635. C           -----------------------
  2636. C
  2637. C           E0     - THE 4 ELEMENTS ON WHICH THE
  2638. C           E1       COMPUTATION OF A NEW ELEMENT IN
  2639. C           E2       THE EPSILON TABLE IS BASED
  2640. C           E3                 E0
  2641. C                        E3    E1    NEW
  2642. C                              E2
  2643. C           NEWELM - NUMBER OF ELEMENTS TO BE COMPUTED IN THE NEW
  2644. C                    DIAGONAL
  2645. C           ERROR  - ERROR = ABS(E1-E0)+ABS(E2-E1)+ABS(NEW-E2)
  2646. C           RESULT - THE ELEMENT IN THE NEW DIAGONAL WITH LEAST VALUE
  2647. C                    OF ERROR
  2648. C
  2649. C           MACHINE DEPENDENT CONSTANTS
  2650. C           ---------------------------
  2651. C
  2652. C           EPMACH IS THE LARGEST RELATIVE SPACING.
  2653. C           OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
  2654. C           LIMEXP IS THE MAXIMUM NUMBER OF ELEMENTS THE EPSILON
  2655. C           TABLE CAN CONTAIN. IF THIS NUMBER IS REACHED, THE UPPER
  2656. C           DIAGONAL OF THE EPSILON TABLE IS DELETED.
  2657. C
  2658. C***FIRST EXECUTABLE STATEMENT  QELG
  2659.       EPMACH = R1MACH(4)
  2660.       OFLOW = R1MACH(2)
  2661.       NRES = NRES+1
  2662.       ABSERR = OFLOW
  2663.       RESULT = EPSTAB(N)
  2664.       IF(N.LT.3) GO TO 100
  2665.       LIMEXP = 50
  2666.       EPSTAB(N+2) = EPSTAB(N)
  2667.       NEWELM = (N-1)/2
  2668.       EPSTAB(N) = OFLOW
  2669.       NUM = N
  2670.       K1 = N
  2671.       DO 40 I = 1,NEWELM
  2672.         K2 = K1-1
  2673.         K3 = K1-2
  2674.         RES = EPSTAB(K1+2)
  2675.         E0 = EPSTAB(K3)
  2676.         E1 = EPSTAB(K2)
  2677.         E2 = RES
  2678.         E1ABS = ABS(E1)
  2679.         DELTA2 = E2-E1
  2680.         ERR2 = ABS(DELTA2)
  2681.         TOL2 = AMAX1(ABS(E2),E1ABS)*EPMACH
  2682.         DELTA3 = E1-E0
  2683.         ERR3 = ABS(DELTA3)
  2684.         TOL3 = AMAX1(E1ABS,ABS(E0))*EPMACH
  2685.         IF(ERR2.GT.TOL2.OR.ERR3.GT.TOL3) GO TO 10
  2686. C
  2687. C           IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE
  2688. C           ACCURACY, CONVERGENCE IS ASSUMED.
  2689. C           RESULT = E2
  2690. C           ABSERR = ABS(E1-E0)+ABS(E2-E1)
  2691. C
  2692.         RESULT = RES
  2693.         ABSERR = ERR2+ERR3
  2694. C ***JUMP OUT OF DO-LOOP
  2695.         GO TO 100
  2696.    10   E3 = EPSTAB(K1)
  2697.         EPSTAB(K1) = E1
  2698.         DELTA1 = E1-E3
  2699.         ERR1 = ABS(DELTA1)
  2700.         TOL1 = AMAX1(E1ABS,ABS(E3))*EPMACH
  2701. C
  2702. C           IF TWO ELEMENTS ARE VERY CLOSE TO EACH OTHER, OMIT
  2703. C           A PART OF THE TABLE BY ADJUSTING THE VALUE OF N
  2704. C
  2705.         IF(ERR1.LE.TOL1.OR.ERR2.LE.TOL2.OR.ERR3.LE.TOL3) GO TO 20
  2706.         SS = 0.1E+01/DELTA1+0.1E+01/DELTA2-0.1E+01/DELTA3
  2707.         EPSINF = ABS(SS*E1)
  2708. C
  2709. C           TEST TO DETECT IRREGULAR BEHAVIOUR IN THE TABLE, AND
  2710. C           EVENTUALLY OMIT A PART OF THE TABLE ADJUSTING THE VALUE
  2711. C           OF N.
  2712. C
  2713.         IF(EPSINF.GT.0.1E-03) GO TO 30
  2714.    20   N = I+I-1
  2715. C ***JUMP OUT OF DO-LOOP
  2716.         GO TO 50
  2717. C
  2718. C           COMPUTE A NEW ELEMENT AND EVENTUALLY ADJUST
  2719. C           THE VALUE OF RESULT.
  2720. C
  2721.    30   RES = E1+0.1E+01/SS
  2722.         EPSTAB(K1) = RES
  2723.         K1 = K1-2
  2724.         ERROR = ERR2+ABS(RES-E2)+ERR3
  2725.         IF(ERROR.GT.ABSERR) GO TO 40
  2726.         ABSERR = ERROR
  2727.         RESULT = RES
  2728.    40 CONTINUE
  2729. C
  2730. C           SHIFT THE TABLE.
  2731. C
  2732.    50 IF(N.EQ.LIMEXP) N = 2*(LIMEXP/2)-1
  2733.       IB = 1
  2734.       IF((NUM/2)*2.EQ.NUM) IB = 2
  2735.       IE = NEWELM+1
  2736.       DO 60 I=1,IE
  2737.         IB2 = IB+2
  2738.         EPSTAB(IB) = EPSTAB(IB2)
  2739.         IB = IB2
  2740.    60 CONTINUE
  2741.       IF(NUM.EQ.N) GO TO 80
  2742.       INDX = NUM-N+1
  2743.       DO 70 I = 1,N
  2744.         EPSTAB(I)= EPSTAB(INDX)
  2745.         INDX = INDX+1
  2746.    70 CONTINUE
  2747.    80 IF(NRES.GE.4) GO TO 90
  2748.       RES3LA(NRES) = RESULT
  2749.       ABSERR = OFLOW
  2750.       GO TO 100
  2751. C
  2752. C           COMPUTE ERROR ESTIMATE
  2753. C
  2754.    90 ABSERR = ABS(RESULT-RES3LA(3))+ABS(RESULT-RES3LA(2))
  2755.      1  +ABS(RESULT-RES3LA(1))
  2756.       RES3LA(1) = RES3LA(2)
  2757.       RES3LA(2) = RES3LA(3)
  2758.       RES3LA(3) = RESULT
  2759.   100 ABSERR = AMAX1(ABSERR,0.5E+01*EPMACH*ABS(RESULT))
  2760.       RETURN
  2761.       END
  2762. *DECK R1MACH
  2763.       REAL FUNCTION R1MACH(I)
  2764. C***BEGIN PROLOGUE  R1MACH
  2765. C***DATE WRITTEN   790101   (YYMMDD)
  2766. C***REVISION DATE  890213   (YYMMDD)
  2767. C***CATEGORY NO.  R1
  2768. C***KEYWORDS  LIBRARY=SLATEC,TYPE=SINGLE PRECISION(R1MACH-S D1MACH-D),
  2769. C             MACHINE CONSTANTS
  2770. C***AUTHOR  FOX, P. A., (BELL LABS)
  2771. C           HALL, A. D., (BELL LABS)
  2772. C           SCHRYER, N. L., (BELL LABS)
  2773. C***PURPOSE  RETURNS SINGLE PRECISION MACHINE DEPENDENT CONSTANTS
  2774. C***DESCRIPTION
  2775. C
  2776. C   R1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS
  2777. C   FOR THE LOCAL MACHINE ENVIRONMENT.  IT IS A FUNCTION
  2778. C   SUBROUTINE WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED
  2779. C   AS FOLLOWS, FOR EXAMPLE
  2780. C
  2781. C        A = R1MACH(I)
  2782. C
  2783. C   WHERE I=1,...,5.  THE (OUTPUT) VALUE OF A ABOVE IS
  2784. C   DETERMINED BY THE (INPUT) VALUE OF I.  THE RESULTS FOR
  2785. C   VARIOUS VALUES OF I ARE DISCUSSED BELOW.
  2786. C
  2787. C   SINGLE-PRECISION MACHINE CONSTANTS
  2788. C   R1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
  2789. C   R1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
  2790. C   R1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING.
  2791. C   R1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING.
  2792. C   R1MACH(5) = LOG10(B)
  2793. C
  2794. C   ASSUME SINGLE PRECISION NUMBERS ARE REPRESENTED IN THE T-DIGIT,
  2795. C   BASE-B FORM
  2796. C
  2797. C              SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
  2798. C
  2799. C   WHERE 0 .LE. X(I) .LT. B FOR I=1,...,T, 0 .LT. X(1), AND
  2800. C   EMIN .LE. E .LE. EMAX.
  2801. C
  2802. C   THE VALUES OF B, T, EMIN AND EMAX ARE PROVIDED IN I1MACH AS
  2803. C   FOLLOWS:
  2804. C   I1MACH(10) = B, THE BASE.
  2805. C   I1MACH(11) = T, THE NUMBER OF BASE-B DIGITS.
  2806. C   I1MACH(12) = EMIN, THE SMALLEST EXPONENT E.
  2807. C   I1MACH(13) = EMAX, THE LARGEST EXPONENT E.
  2808. C
  2809. C   TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT,
  2810. C   THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY
  2811. C   REMOVING THE C FROM COLUMN 1.  ALSO, THE VALUES OF
  2812. C   R1MACH(1) - R1MACH(4) SHOULD BE CHECKED FOR CONSISTENCY
  2813. C   WITH THE LOCAL OPERATING SYSTEM.
  2814. C
  2815. C***REFERENCES  FOX, P.A., HALL, A.D., SCHRYER, N.L, *FRAMEWORK FOR
  2816. C                 A PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHE-
  2817. C                 MATICAL SOFTWARE, VOL. 4, NO. 2, JUNE 1978,
  2818. C                 PP. 177-188.
  2819. C***ROUTINES CALLED  XERROR
  2820. C
  2821. C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
  2822. C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
  2823. C        NOTE ADDED BY F. ROMANI  7/11/89
  2824. C
  2825. C***END PROLOGUE  R1MACH
  2826. C
  2827.       INTEGER SMALL(2)
  2828.       INTEGER LARGE(2)
  2829.       INTEGER RIGHT(2)
  2830.       INTEGER DIVER(2)
  2831.       INTEGER LOG10(2)
  2832. C
  2833.       REAL RMACH(5)
  2834.       SAVE RMACH
  2835. C
  2836.       EQUIVALENCE (RMACH(1),SMALL(1))
  2837.       EQUIVALENCE (RMACH(2),LARGE(1))
  2838.       EQUIVALENCE (RMACH(3),RIGHT(1))
  2839.       EQUIVALENCE (RMACH(4),DIVER(1))
  2840.       EQUIVALENCE (RMACH(5),LOG10(1))
  2841. C
  2842. C     MACHINE CONSTANTS FOR THE AMIGA
  2843. C     ABSOFT FORTRAN COMPILER USING THE 68020/68881 COMPILER OPTION
  2844. C
  2845. C     DATA SMALL(1) / Z'00800000' /
  2846. C     DATA LARGE(1) / Z'7F7FFFFF' /
  2847. C     DATA RIGHT(1) / Z'33800000' /
  2848. C     DATA DIVER(1) / Z'34000000' /
  2849. C     DATA LOG10(1) / Z'3E9A209B' /
  2850. C
  2851. C     MACHINE CONSTANTS FOR THE AMIGA
  2852. C     ABSOFT FORTRAN COMPILER USING SOFTWARE FLOATING POINT
  2853. C
  2854. C     DATA SMALL(1) / Z'00800000' /
  2855. C     DATA LARGE(1) / Z'7EFFFFFF' /
  2856. C     DATA RIGHT(1) / Z'33800000' /
  2857. C     DATA DIVER(1) / Z'34000000' /
  2858. C     DATA LOG10(1) / Z'3E9A209B' /
  2859. C
  2860. C     MACHINE CONSTANTS FOR THE APOLLO
  2861. C
  2862. C     DATA SMALL(1) / 16#00800000 /
  2863. C     DATA LARGE(1) / 16#7FFFFFFF /
  2864. C     DATA RIGHT(1) / 16#33800000 /
  2865. C     DATA DIVER(1) / 16#34000000 /
  2866. C     DATA LOG10(1) / 16#3E9A209B /
  2867. C
  2868. C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM
  2869. C
  2870. C     DATA RMACH(1) / Z400800000 /
  2871. C     DATA RMACH(2) / Z5FFFFFFFF /
  2872. C     DATA RMACH(3) / Z4E9800000 /
  2873. C     DATA RMACH(4) / Z4EA800000 /
  2874. C     DATA RMACH(5) / Z500E730E8 /
  2875. C
  2876. C     MACHINE CONSTANTS FOR THE BURROUGHS 5700/6700/7700 SYSTEMS
  2877. C
  2878. C     DATA RMACH(1) / O1771000000000000 /
  2879. C     DATA RMACH(2) / O0777777777777777 /
  2880. C     DATA RMACH(3) / O1311000000000000 /
  2881. C     DATA RMACH(4) / O1301000000000000 /
  2882. C     DATA RMACH(5) / O1157163034761675 /
  2883. C
  2884. C     MACHINE CONSTANTS FOR THE CDC 170/180 SERIES USING NOS/VE
  2885. C
  2886. C     DATA RMACH(1) / Z"3001800000000000" /
  2887. C     DATA RMACH(2) / Z"4FFEFFFFFFFFFFFE" /
  2888. C     DATA RMACH(3) / Z"3FD2800000000000" /
  2889. C     DATA RMACH(4) / Z"3FD3800000000000" /
  2890. C     DATA RMACH(5) / Z"3FFF9A209A84FBCF" /
  2891. C
  2892. C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES
  2893. C
  2894. C     DATA RMACH(1) / 00564000000000000000B /
  2895. C     DATA RMACH(2) / 37767777777777777776B /
  2896. C     DATA RMACH(3) / 16414000000000000000B /
  2897. C     DATA RMACH(4) / 16424000000000000000B /
  2898. C     DATA RMACH(5) / 17164642023241175720B /
  2899. C
  2900. C     MACHINE CONSTANTS FOR THE CELERITY C1260
  2901. C
  2902. C     DATA SMALL(1) / Z'00800000' /
  2903. C     DATA LARGE(1) / Z'7F7FFFFF' /
  2904. C     DATA RIGHT(1) / Z'33800000' /
  2905. C     DATA DIVER(1) / Z'34000000' /
  2906. C     DATA LOG10(1) / Z'3E9A209B' /
  2907. C
  2908. C     MACHINE CONSTANTS FOR THE CONVEX C-1
  2909. C
  2910. C     DATA SMALL(1) / '00800000'X /
  2911. C     DATA LARGE(1) / '7FFFFFFF'X /
  2912. C     DATA RIGHT(1) / '34800000'X /
  2913. C     DATA DIVER(1) / '35000000'X /
  2914. C     DATA LOG10(1) / '3F9A209B'X /
  2915. C
  2916. C     MACHINE CONSTANTS FOR THE CRAY-1
  2917. C
  2918. C     DATA RMACH(1) / 200034000000000000000B /
  2919. C     DATA RMACH(2) / 577767777777777777776B /
  2920. C     DATA RMACH(3) / 377224000000000000000B /
  2921. C     DATA RMACH(4) / 377234000000000000000B /
  2922. C     DATA RMACH(5) / 377774642023241175720B /
  2923. C
  2924. C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
  2925. C
  2926. C     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD -
  2927. C     STATIC RMACH(5)
  2928. C
  2929. C     DATA SMALL /    20K,       0 /
  2930. C     DATA LARGE / 77777K, 177777K /
  2931. C     DATA RIGHT / 35420K,       0 /
  2932. C     DATA DIVER / 36020K,       0 /
  2933. C     DATA LOG10 / 40423K,  42023K /
  2934. C
  2935. C     MACHINE CONSTANTS FOR THE HARRIS 220
  2936. C
  2937. C     DATA SMALL(1), SMALL(2) / '20000000, '00000201 /
  2938. C     DATA LARGE(1), LARGE(2) / '37777777, '00000177 /
  2939. C     DATA RIGHT(1), RIGHT(2) / '20000000, '00000352 /
  2940. C     DATA DIVER(1), DIVER(2) / '20000000, '00000353 /
  2941. C     DATA LOG10(1), LOG10(2) / '23210115, '00000377 /
  2942. C
  2943. C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES
  2944. C
  2945. C     DATA RMACH(1) / O402400000000 /
  2946. C     DATA RMACH(2) / O376777777777 /
  2947. C     DATA RMACH(3) / O714400000000 /
  2948. C     DATA RMACH(4) / O716400000000 /
  2949. C     DATA RMACH(5) / O776464202324 /
  2950. C
  2951. C     MACHINE CONSTANTS FOR THE HP 2100
  2952. C     3 WORD DOUBLE PRECISION WITH FTN4
  2953. C
  2954. C     DATA SMALL(1), SMALL(2) / 40000B,       1 /
  2955. C     DATA LARGE(1), LARGE(2) / 77777B, 177776B /
  2956. C     DATA RIGHT(1), RIGHT(2) / 40000B,    325B /
  2957. C     DATA DIVER(1), DIVER(2) / 40000B,    327B /
  2958. C     DATA LOG10(1), LOG10(2) / 46420B,  46777B /
  2959. C
  2960. C     MACHINE CONSTANTS FOR THE HP 2100
  2961. C     4 WORD DOUBLE PRECISION WITH FTN4
  2962. C
  2963. C     DATA SMALL(1), SMALL(2) / 40000B,       1 /
  2964. C     DATA LARGE91), LARGE(2) / 77777B, 177776B /
  2965. C     DATA RIGHT(1), RIGHT(2) / 40000B,    325B /
  2966. C     DATA DIVER(1), DIVER(2) / 40000B,    327B /
  2967. C     DATA LOG10(1), LOG10(2) / 46420B,  46777B /
  2968. C
  2969. C     MACHINE CONSTANTS FOR THE HP 9000
  2970. C
  2971. C     DATA SMALL(1) / 00004000000B /
  2972. C     DATA LARGE(1) / 17677777777B /
  2973. C     DATA RIGHT(1) / 06340000000B /
  2974. C     DATA DIVER(1) / 06400000000B /
  2975. C     DATA LOG10(1) / 07646420233B /
  2976. C
  2977. C     MACHINE CONSTANTS FOR THE ELXSI 6400
  2978. C       ASSUMING REAL*4 IS THE DEFAULT REAL
  2979. C
  2980. C     DATA SMALL(1) / '00800000'X /
  2981. C     DATA LARGE(1) / '7F7FFFFF'X /
  2982. C     DATA RIGHT(1) / '33800000'X /
  2983. C     DATA DIVER(1) / '34000000'X /
  2984. C     DATA LOG10(1) / '3E9A209B'X /
  2985. C
  2986. C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
  2987. C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86  AND
  2988. C     THE PERKIN ELMER (INTERDATA) 7/32.
  2989. C
  2990. C     DATA RMACH(1) / Z00100000 /
  2991. C     DATA RMACH(2) / Z7FFFFFFF /
  2992. C     DATA RMACH(3) / Z3B100000 /
  2993. C     DATA RMACH(4) / Z3C100000 /
  2994. C     DATA RMACH(5) / Z41134413 /
  2995. C
  2996. C     MACHINE CONSTANTS FOR THE IBM PC
  2997. C
  2998. C     DATA SMALL(1) / 1.18E-38      /
  2999. C     DATA LARGE(1) / 3.40E+38      /
  3000. C     DATA RIGHT(1) / 0.595E-07     /
  3001. C     DATA DIVER(1) / 1.19E-07      /
  3002. C     DATA LOG10(1) / 0.30102999566 /
  3003. C
  3004. C     MACHINE CONSTANTS FOR THE PDP-10 (KA OR KI PROCESSOR)
  3005. C
  3006. C     DATA RMACH(1) / "000400000000 /
  3007. C     DATA RMACH(2) / "377777777777 /
  3008. C     DATA RMACH(3) / "146400000000 /
  3009. C     DATA RMACH(4) / "147400000000 /
  3010. C     DATA RMACH(5) / "177464202324 /
  3011. C
  3012. C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
  3013. C     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
  3014. C
  3015. C     DATA SMALL(1) /    8388608 /
  3016. C     DATA LARGE(1) / 2147483647 /
  3017. C     DATA RIGHT(1) /  880803840 /
  3018. C     DATA DIVER(1) /  889192448 /
  3019. C     DATA LOG10(1) / 1067065499 /
  3020. C
  3021. C     DATA RMACH(1) / O00040000000 /
  3022. C     DATA RMACH(2) / O17777777777 /
  3023. C     DATA RMACH(3) / O06440000000 /
  3024. C     DATA RMACH(4) / O06500000000 /
  3025. C     DATA RMACH(5) / O07746420233 /
  3026. C
  3027. C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
  3028. C     16-BIT INTEGERS  (EXPRESSED IN INTEGER AND OCTAL).
  3029. C
  3030. C     DATA SMALL(1), SMALL(2) /   128,     0 /
  3031. C     DATA LARGE(1), LARGE(2) / 32767,    -1 /
  3032. C     DATA RIGHT(1), RIGHT(2) / 13440,     0 /
  3033. C     DATA DIVER(1), DIVER(2) / 13568,     0 /
  3034. C     DATA LOG10(1), LOG10(2) / 16282,  8347 /
  3035. C
  3036. C     DATA SMALL(1), SMALL(2) / O000200, O000000 /
  3037. C     DATA LARGE(1), LARGE(2) / O077777, O177777 /
  3038. C     DATA RIGHT(1), RIGHT(2) / O032200, O000000 /
  3039. C     DATA DIVER(1), DIVER(2) / O032400, O000000 /
  3040. C     DATA LOG10(1), LOG10(2) / O037632, O020233 /
  3041. C
  3042. C     MACHINE CONSTANTS FOR THE SUN
  3043. C
  3044. C     DATA SMALL(1) / Z'00800000' /
  3045. C     DATA LARGE(1) / Z'7F7FFFFF' /
  3046. C     DATA RIGHT(1) / Z'33800000' /
  3047. C     DATA DIVER(1) / Z'34000000' /
  3048. C     DATA LOG10(1) / Z'3E9A209B' /
  3049. C
  3050. C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES
  3051. C
  3052. C     DATA RMACH(1) / O000400000000 /
  3053. C     DATA RMACH(2) / O377777777777 /
  3054. C     DATA RMACH(3) / O146400000000 /
  3055. C     DATA RMACH(4) / O147400000000 /
  3056. C     DATA RMACH(5) / O177464202324 /
  3057. C
  3058. C     MACHINE CONSTANTS FOR THE VAX 11/780
  3059. C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
  3060. C     THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS
  3061. C     THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS
  3062. C
  3063. C     DATA SMALL(1) /       128 /
  3064. C     DATA LARGE(1) /    -32769 /
  3065. C     DATA RIGHT(1) /     13440 /
  3066. C     DATA DIVER(1) /     13568 /
  3067. C     DATA LOG10(1) / 547045274 /
  3068. C
  3069. C     DATA SMALL(1) / Z00000080 /
  3070. C     DATA LARGE(1) / ZFFFF7FFF /
  3071. C     DATA RIGHT(1) / Z00003480 /
  3072. C     DATA DIVER(1) / Z00003500 /
  3073. C     DATA LOG10(1) / Z209B3F9A /
  3074. C
  3075. C     MACHINE CONSTANTS FOR THE Z80 MICROPROCESSOR
  3076. C
  3077. C     DATA SMALL(1), SMALL(2) /     0,    256/
  3078. C     DATA LARGE(1), LARGE(2) /    -1,   -129/
  3079. C     DATA RIGHT(1), RIGHT(2) /     0,  26880/
  3080. C     DATA DIVER(1), DIVER(2) /     0,  27136/
  3081. C     DATA LOG10(1), LOG10(2) /  8347,  32538/
  3082. C
  3083. C
  3084. C***FIRST EXECUTABLE STATEMENT  R1MACH
  3085. C
  3086. C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
  3087. C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
  3088. C        NOTE ADDED BY F. ROMANI  7/11/89
  3089. C
  3090. C     IF (I .LT. 1  .OR.  I .GT. 5)
  3091. C    1   CALL XERROR ('R1MACH -- I OUT OF BOUNDS', 25, 1, 2)
  3092. C
  3093.       R1MACH = RMACH(I)
  3094.       RETURN
  3095. C
  3096.       END
  3097. *DECK QXGS
  3098.       SUBROUTINE QXGS(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,IER,
  3099.      1   LIMIT,LENIW,LENW,LAST,IWORK,WORK)
  3100. C
  3101. C            THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
  3102. C            DEFINITE INTEGRAL  I = INTEGRAL OF F OVER (A,B),
  3103. C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
  3104. C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
  3105. C
  3106. C        PARAMETERS
  3107. C         ON ENTRY
  3108. C            F      - REAL
  3109. C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  3110. C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  3111. C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  3112. C
  3113. C            A      - REAL
  3114. C                     LOWER LIMIT OF INTEGRATION
  3115. C
  3116. C            B      - REAL
  3117. C                     UPPER LIMIT OF INTEGRATION
  3118. C
  3119. C            EPSABS - REAL
  3120. C                     ABSOLUTE ACCURACY REQUESTED
  3121. C            EPSREL - REAL
  3122. C                     RELATIVE ACCURACY REQUESTED
  3123. C                     IF  EPSABS.LE.0
  3124. C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28),
  3125. C                     THE ROUTINE WILL END WITH IER = 6.
  3126. C
  3127. C         ON RETURN
  3128. C            RESULT - REAL
  3129. C                     APPROXIMATION TO THE INTEGRAL
  3130. C
  3131. C            ABSERR - REAL
  3132. C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
  3133. C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
  3134. C
  3135. C            IER    - INTEGER
  3136. C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
  3137. C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
  3138. C                             ACCURACY HAS BEEN ACHIEVED.
  3139. C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE
  3140. C                             THE ESTIMATES FOR INTEGRAL AND ERROR ARE
  3141. C                             LESS RELIABLE. IT IS ASSUMED THAT THE
  3142. C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
  3143. C            ERROR MESSAGES
  3144. C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
  3145. C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB-
  3146. C                             DIVISIONS BY INCREASING THE VALUE OF LIMIT
  3147. C                             (AND TAKING THE ACCORDING DIMENSION
  3148. C                             ADJUSTMENTS INTO ACCOUNT. HOWEVER, IF
  3149. C                             THIS YIELDS NO IMPROVEMENT IT IS ADVISED
  3150. C                             TO ANALYZE THE INTEGRAND IN ORDER TO
  3151. C                             DETERMINE THE INTEGRATION DIFFICULTIES. IF
  3152. C                             THE POSITION OF A LOCAL DIFFICULTY CAN BE
  3153. C                             DETERMINED (E.G. SINGULARITY,
  3154. C                             DISCONTINUITY WITHIN THE INTERVAL) ONE
  3155. C                             WILL PROBABLY GAIN FROM SPLITTING UP THE
  3156. C                             INTERVAL AT THIS POINT AND CALLING THE
  3157. C                             INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
  3158. C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
  3159. C                             SHOULD BE USED, WHICH IS DESIGNED FOR
  3160. C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.
  3161. C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS DETEC-
  3162. C                             TED, WHICH PREVENTS THE REQUESTED
  3163. C                             TOLERANCE FROM BEING ACHIEVED.
  3164. C                             THE ERROR MAY BE UNDER-ESTIMATED.
  3165. C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR
  3166. C                             OCCURS AT SOME POINTS OF THE INTEGRATION
  3167. C                             INTERVAL.
  3168. C                         = 4 THE ALGORITHM DOES NOT CONVERGE.
  3169. C                             ROUNDOFF ERROR IS DETECTED IN THE
  3170. C                             EXTRAPOLATION TABLE. IT IS PRESUMED THAT
  3171. C                             THE REQUESTED TOLERANCE CANNOT BE
  3172. C                             ACHIEVED, AND THAT THE RETURNED RESULT IS
  3173. C                             THE BEST WHICH CAN BE OBTAINED.
  3174. C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
  3175. C                             SLOWLY CONVERGENT. IT MUST BE NOTED THAT
  3176. C                             DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
  3177. C                             OF IER.
  3178. C                         = 6 THE INPUT IS INVALID, BECAUSE
  3179. C                             (EPSABS.LE.0 AND
  3180. C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28)
  3181. C                             OR LIMIT.LT.1 OR LENW.LT.LIMIT*46 OR
  3182. C                                LENIW .LT. LIMIT*3
  3183. C                             RESULT, ABSERR, LAST ARE SET TO
  3184. C                             ZERO.EXCEPT WHEN LIMIT OR LENW OR LENIW
  3185. C                             IS INVALID, IWORK(1), WORK(LIMIT*2+1) AND
  3186. C                             WORK(LIMIT*3+1) ARE SET TO ZERO, WORK(1)
  3187. C                             IS SET TO A AND WORK(LIMIT+1) TO B.
  3188. C
  3189. C         DIMENSIONING PARAMETERS
  3190. C            LIMIT - INTEGER
  3191. C                    LIMIT DETERMINES THE MAXIMUM NUMBER OF SUBINTERVALS
  3192. C                    IN THE PARTITION OF THE GIVEN INTEGRATION INTERVAL
  3193. C                    (A,B), LIMIT.GE.1.
  3194. C                    IF LIMIT.LT.1, THE ROUTINE WILL END WITH IER = 6.
  3195. C
  3196. C            LENW  - INTEGER
  3197. C                    DIMENSIONING PARAMETER FOR WORK
  3198. C                    LENW MUST BE AT LEAST LIMIT*46
  3199. C                    IF LENW.LT.LIMIT*46, THE ROUTINE WILL END
  3200. C                    WITH IER = 6.
  3201. C
  3202. C            LENIW - INTEGER
  3203. C                    DIMENSIONING PARAMETER FOR IWORK
  3204. C                    LENIW MUST BE AT LEAST LIMIT*3
  3205. C                    IF LENW.LT.LIMIT*3, THE ROUTINE WILL END
  3206. C                    WITH IER = 6.
  3207. C
  3208. C            LAST  - INTEGER
  3209. C                    ON RETURN, LAST EQUALS THE NUMBER OF SUBINTERVALS
  3210. C                    PRODUCED IN THE SUBDIVISION PROCESS, DETERMINES THE
  3211. C                    NUMBER OF SIGNIFICANT ELEMENTS ACTUALLY IN THE WORK
  3212. C                    ARRAYS.
  3213. C
  3214. C         WORK ARRAYS
  3215. C            IWORK - INTEGER
  3216. C                    VECTOR OF DIMENSION AT LEAST 3*LIMIT, THE FIRST K
  3217. C                    ELEMENTS OF WHICH CONTAIN POINTERS
  3218. C                    TO THE ERROR ESTIMATES OVER THE SUBINTERVALS
  3219. C                    SUCH THAT WORK(LIMIT*3+IWORK(1)),... ,
  3220. C                    WORK(LIMIT*3+IWORK(K)) FORM A DECREASING
  3221. C                    SEQUENCE, WITH K = LAST IF LAST.LE.(LIMIT/2+2),
  3222. C                    AND K = LIMIT+1-LAST OTHERWISE
  3223. C
  3224. C            WORK  - REAL
  3225. C                    VECTOR OF DIMENSION AT LEAST LENW
  3226. C                    ON RETURN
  3227. C                    WORK(1), ..., WORK(LAST) CONTAIN THE LEFT
  3228. C                     END-POINTS OF THE SUBINTERVALS IN THE
  3229. C                     PARTITION OF (A,B),
  3230. C                    WORK(LIMIT+1), ..., WORK(LIMIT+LAST) CONTAIN
  3231. C                     THE RIGHT END-POINTS,
  3232. C                    WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) CONTAIN
  3233. C                     THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS,
  3234. C                    WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST)
  3235. C                     CONTAIN THE ERROR ESTIMATES.
  3236. C                    WORK(LIMIT*4+1), ... IS THE AREA RESERVED TO STORE
  3237. C                     FUNCTIONAL VALUES.
  3238. C***ROUTINES CALLED  QXGSE,XERROR
  3239. C
  3240. C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
  3241. C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
  3242. C
  3243.       REAL A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK
  3244.       INTEGER IER,IWORK,LAST,LENIW,LENW,LIMIT,LVL,L1,L2,L3,L4,L5
  3245. C
  3246.       DIMENSION IWORK(LENIW),WORK(LENW)
  3247. C
  3248.       EXTERNAL F
  3249. C
  3250. C         CHECK VALIDITY OF LIMIT,LENIW AND LENW.
  3251. C
  3252. C***FIRST EXECUTABLE STATEMENT  QXGS
  3253.       IER = 6
  3254.       LAST = 0
  3255.       RESULT = 0.0E+00
  3256.       ABSERR = 0.0E+00
  3257.       IF(LIMIT.LT.1.OR.LENIW.LT.LIMIT*3.OR.LENW.LT.LIMIT*46) GO TO 10
  3258. C
  3259. C         PREPARE CALL FOR QXGSE.
  3260. C
  3261.       L1 = LIMIT+1
  3262.       L2 = LIMIT+L1
  3263.       L3 = LIMIT+L2
  3264.       L4 = LIMIT+L3
  3265.       L5 = 21*LIMIT+L3
  3266. C
  3267.       CALL QXGSE(F,A,B,EPSABS,EPSREL,LIMIT,RESULT,ABSERR,
  3268.      *  IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST,
  3269.      *  WORK(L4),WORK(L5),IWORK(L1),IWORK(L2))
  3270. C
  3271. C         CALL ERROR HANDLER IF NECESSARY.
  3272. C
  3273.       LVL = 0
  3274. 10    IF(IER.EQ.6) LVL = 1
  3275. C
  3276. C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
  3277. C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
  3278. C
  3279. C     IF(IER.NE.0)CALL XERROR('ABNORMAL RETURN FROM QXGS  ',28,IER,LVL)
  3280.       RETURN
  3281.       END
  3282. *DECK QXGSE
  3283.       SUBROUTINE QXGSE(F,A,B,EPSABS,EPSREL,LIMIT,RESULT,ABSERR,
  3284.      *   IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST,VALP,VALN,LP,LN)
  3285. C
  3286. C            THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
  3287. C            DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B),
  3288. C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
  3289. C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
  3290. C
  3291. C        PARAMETERS
  3292. C         ON ENTRY
  3293. C            F      - REAL
  3294. C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  3295. C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  3296. C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  3297. C
  3298. C            A      - REAL
  3299. C                     LOWER LIMIT OF INTEGRATION
  3300. C
  3301. C            B      - REAL
  3302. C                     UPPER LIMIT OF INTEGRATION
  3303. C
  3304. C            EPSABS - REAL
  3305. C                     ABSOLUTE ACCURACY REQUESTED
  3306. C            EPSREL - REAL
  3307. C                     RELATIVE ACCURACY REQUESTED
  3308. C                     IF  EPSABS.LE.0
  3309. C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28),
  3310. C                     THE ROUTINE WILL END WITH IER = 6.
  3311. C
  3312. C            LIMIT  - INTEGER
  3313. C                     GIVES AN UPPERBOUND ON THE NUMBER OF SUBINTERVALS
  3314. C                     IN THE PARTITION OF (A,B)
  3315. C
  3316. C         ON RETURN
  3317. C            RESULT - REAL
  3318. C                     APPROXIMATION TO THE INTEGRAL
  3319. C
  3320. C            ABSERR - REAL
  3321. C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
  3322. C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
  3323. C
  3324. C            IER    - INTEGER
  3325. C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
  3326. C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
  3327. C                             ACCURACY HAS BEEN ACHIEVED.
  3328. C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE
  3329. C                             THE ESTIMATES FOR INTEGRAL AND ERROR ARE
  3330. C                             LESS RELIABLE. IT IS ASSUMED THAT THE
  3331. C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
  3332. C            ERROR MESSAGES
  3333. C                         = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
  3334. C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB-
  3335. C                             DIVISIONS BY INCREASING THE VALUE OF LIMIT
  3336. C                             (AND TAKING THE ACCORDING DIMENSION
  3337. C                             ADJUSTMENTS INTO ACCOUNT). HOWEVER, IF
  3338. C                             THIS YIELDS NO IMPROVEMENT IT IS ADVISED
  3339. C                             TO ANALYZE THE INTEGRAND IN ORDER TO
  3340. C                             DETERMINE THE INTEGRATION DIFFICULTIES. IF
  3341. C                             THE POSITION OF A LOCAL DIFFICULTY CAN BE
  3342. C                             DETERMINED (E.G. SINGULARITY,
  3343. C                             DISCONTINUITY WITHIN THE INTERVAL) ONE
  3344. C                             WILL PROBABLY GAIN FROM SPLITTING UP THE
  3345. C                             INTERVAL AT THIS POINT AND CALLING THE
  3346. C                             INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
  3347. C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
  3348. C                             SHOULD BE USED, WHICH IS DESIGNED FOR
  3349. C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.
  3350. C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS DETEC-
  3351. C                             TED, WHICH PREVENTS THE REQUESTED
  3352. C                             TOLERANCE FROM BEING ACHIEVED.
  3353. C                             THE ERROR MAY BE UNDER-ESTIMATED.
  3354. C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR
  3355. C                             OCCURS AT SOME POINTS OF THE INTEGRATION
  3356. C                             INTERVAL.
  3357. C                         = 4 THE ALGORITHM DOES NOT CONVERGE.
  3358. C                             ROUNDOFF ERROR IS DETECTED IN THE
  3359. C                             EXTRAPOLATION TABLE.
  3360. C                             IT IS PRESUMED THAT THE REQUESTED
  3361. C                             TOLERANCE CANNOT BE ACHIEVED, AND THAT THE
  3362. C                             RETURNED RESULT IS THE BEST WHICH CAN BE
  3363. C                             OBTAINED.
  3364. C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
  3365. C                             SLOWLY CONVERGENT. IT MUST BE NOTED THAT
  3366. C                             DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
  3367. C                             OF IER.
  3368. C                         = 6 THE INPUT IS INVALID, BECAUSE
  3369. C                             EPSABS.LE.0 AND
  3370. C                             EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28).
  3371. C                             RESULT, ABSERR, LAST, RLIST(1),
  3372. C                             IORD(1) AND ELIST(1) ARE SET TO ZERO.
  3373. C                             ALIST(1) AND BLIST(1) ARE SET TO A AND B
  3374. C                             RESPECTIVELY.
  3375. C
  3376. C            ALIST  - REAL
  3377. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  3378. C                      LAST  ELEMENTS OF WHICH ARE THE LEFT END POINTS
  3379. C                     OF THE SUBINTERVALS IN THE PARTITION OF THE
  3380. C                     GIVEN INTEGRATION RANGE (A,B)
  3381. C
  3382. C            BLIST  - REAL
  3383. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  3384. C                      LAST  ELEMENTS OF WHICH ARE THE RIGHT END POINTS
  3385. C                     OF THE SUBINTERVALS IN THE PARTITION OF THE GIVEN
  3386. C                     INTEGRATION RANGE (A,B)
  3387. C
  3388. C            RLIST  - REAL
  3389. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  3390. C                      LAST  ELEMENTS OF WHICH ARE THE INTEGRAL
  3391. C                     APPROXIMATIONS ON THE SUBINTERVALS
  3392. C
  3393. C            ELIST  - REAL
  3394. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  3395. C                      LAST  ELEMENTS OF WHICH ARE THE MODULI OF THE
  3396. C                     ABSOLUTE ERROR ESTIMATES ON THE SUBINTERVALS
  3397. C
  3398. C            IORD   - INTEGER
  3399. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST K
  3400. C                     ELEMENTS OF WHICH ARE POINTERS TO THE
  3401. C                     ERROR ESTIMATES OVER THE SUBINTERVALS,
  3402. C                     SUCH THAT ELIST(IORD(1)), ..., ELIST(IORD(K))
  3403. C                     FORM A DECREASING SEQUENCE, WITH K = LAST
  3404. C                     IF LAST.LE.(LIMIT/2+2), AND K = LIMIT+1-LAST
  3405. C                     OTHERWISE
  3406. C
  3407. C            LAST   - INTEGER
  3408. C                     NUMBER OF SUBINTERVALS ACTUALLY PRODUCED IN THE
  3409. C                     SUBDIVISION PROCESS
  3410. C
  3411. C            VALP   - REAL
  3412. C            VALN     ARRAYS OF DIMENSION AT LEAST (21,LIMIT) USED TO
  3413. C                     SAVE THE FUNCTIONAL VALUES
  3414. C
  3415. C            LP     - INTEGER
  3416. C            LN       VECTORS OF DIMENSION AT LEAST LIMIT, USED TO
  3417. C                     STORE THE ACTUAL NUMBER OF FUNCTIONAL VALUES
  3418. C                     SAVED IN THE CORRESPONDING COLUMN
  3419. C                     OF VALP,VALN
  3420. C
  3421. C***ROUTINES CALLED  F,R1MACH,QELG,QXLQM,QPSRT,QXRRD,QXCPY
  3422. C
  3423.       REAL A,ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,
  3424.      *  A2,B,BLIST,B1,B2,CORREC,ABS,DEFABS,DEFAB1,DEFAB2,R1MACH,AMAX1,
  3425.      *  DRES,ELIST,EPMACH,EPSABS,EPSREL,ERLARG,ERLAST,ERRBND,ERRMAX,
  3426.      *  ERROR1,ERROR2,ERRO12,ERRSUM,ERTEST,F,OFLOW,RESABS,RESEPS,RESULT,
  3427.      *  RES3LA,RLIST,RLIST2,SMALL,UFLOW,
  3428.      *  VALP,VALN,VP1,VP2,VN1,VN2
  3429.       INTEGER ID,IER,IERRO,IORD,IROFF1,IROFF2,IROFF3,JUPBND,K,KSGN,
  3430.      *  KTMIN,LAST,LIMIT,MAXERR,NRES,NRMAX,NUMRL2,
  3431.      *  LP,LN,LP1,LP2,LN1,LN2
  3432.       LOGICAL EXTRAP,NOEXT
  3433.       DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT),
  3434.      * RES3LA(3),RLIST(LIMIT),RLIST2(52),
  3435.      * VALP(21,LIMIT),VALN(21,LIMIT),LP(LIMIT),LN(LIMIT),
  3436.      * VP1(21),VP2(21),VN1(21),VN2(21)
  3437.       EXTERNAL F
  3438. C
  3439. C            MACHINE DEPENDENT CONSTANTS
  3440. C            ---------------------------
  3441. C
  3442. C          EPMACH IS THE LARGEST RELATIVE SPACING.
  3443. C          UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
  3444. C          OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
  3445. C
  3446. C***FIRST EXECUTABLE STATEMENT  QXGSE
  3447.       EPMACH = R1MACH(4)
  3448. C
  3449. C            TEST ON VALIDITY OF PARAMETERS
  3450. C            ------------------------------
  3451.       IER = 0
  3452.       LAST = 0
  3453.       RESULT = 0.0E+00
  3454.       ABSERR = 0.0E+00
  3455.       ALIST(1) = A
  3456.       BLIST(1) = B
  3457.       RLIST(1) = 0.0E+00
  3458.       ELIST(1) = 0.0E+00
  3459.       IF(EPSABS.LE.0.0E+00.AND.EPSREL.LT.AMAX1(0.5E+02*EPMACH,0.5E-28))
  3460.      *   IER = 6
  3461.       IF(IER.EQ.6) GO TO 999
  3462. C
  3463. C           FIRST APPROXIMATION TO THE INTEGRAL
  3464. C           -----------------------------------
  3465. C
  3466.       UFLOW = R1MACH(1)
  3467.       OFLOW = R1MACH(2)
  3468.       IERRO = 0
  3469.       LP(1)=1
  3470.       LN(1)=1
  3471.       VALP(1,1)=F((A+B)*0.5E0)
  3472.       VALN(1,1)=VALP(1,1)
  3473.       CALL QXLQM(F,A,B,RESULT,ABSERR,DEFABS,RESABS,
  3474.      *           VALP(1,1),VALN(1,1),LP(1),LN(1), 2  )
  3475. C
  3476. C           TEST ON ACCURACY.
  3477. C
  3478.       DRES = ABS(RESULT)
  3479.       ERRBND = AMAX1(EPSABS,EPSREL*DRES)
  3480.       LAST = 1
  3481.       RLIST(1) = RESULT
  3482.       ELIST(1) = ABSERR
  3483.       IORD(1) = 1
  3484.       IF(ABSERR.LE.1.0E+02*EPMACH*DEFABS.AND.ABSERR.GT.ERRBND) IER = 2
  3485.       IF(LIMIT.EQ.1) IER = 1
  3486.       IF(IER.NE.0.OR.(ABSERR.LE.ERRBND.AND.ABSERR.NE.RESABS).OR.
  3487.      *  ABSERR.EQ.0.0E+00) GO TO 999
  3488. C
  3489. C           INITIALIZATION
  3490. C           --------------
  3491. C
  3492.       RLIST2(1) = RESULT
  3493.       ERRMAX = ABSERR
  3494.       MAXERR = 1
  3495.       AREA = RESULT
  3496.       ERRSUM = ABSERR
  3497.       ABSERR = OFLOW
  3498.       NRMAX = 1
  3499.       NRES = 0
  3500.       NUMRL2 = 2
  3501.       KTMIN = 0
  3502.       EXTRAP = .FALSE.
  3503.       NOEXT = .FALSE.
  3504.       IROFF1 = 0
  3505.       IROFF2 = 0
  3506.       IROFF3 = 0
  3507.       KSGN = -1
  3508.       IF(DRES.GE.(0.1E+01-0.5E+02*EPMACH)*DEFABS) KSGN = 1
  3509. C
  3510. C           MAIN DO-LOOP
  3511. C           ------------
  3512. C
  3513.       DO 90 LAST = 2,LIMIT
  3514. C
  3515. C           BISECT THE SUBINTERVAL WITH THE NRMAX-TH LARGEST ERROR
  3516. C           ESTIMATE.
  3517. C
  3518.         A1 = ALIST(MAXERR)
  3519.         B1 = 0.5E+00*(ALIST(MAXERR)+BLIST(MAXERR))
  3520.         A2 = B1
  3521.         B2 = BLIST(MAXERR)
  3522.         ERLAST = ERRMAX
  3523.         CALL QXRRD(F,VALN(1,MAXERR),LN(MAXERR),B1,A1,VN1,VP1,LN1,LP1)
  3524.         CALL QXLQM(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1,VP1,VN1,LP1,LN1,
  3525.      *              2)
  3526.         CALL QXRRD(F,VALP(1,MAXERR),LP(MAXERR),A2,B2,VP2,VN2,LP2,LN2)
  3527.         CALL QXLQM(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2,VP2,VN2,LP2,LN2,
  3528.      *              2)
  3529. C
  3530. C           IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
  3531. C           AND ERROR AND TEST FOR ACCURACY.
  3532. C
  3533.         AREA12 = AREA1+AREA2
  3534.         ERRO12 = ERROR1+ERROR2
  3535.         ERRSUM = ERRSUM+ERRO12-ERRMAX
  3536.         AREA = AREA+AREA12-RLIST(MAXERR)
  3537.         IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2) GO TO 15
  3538.         IF(ABS(RLIST(MAXERR)-AREA12).GT.0.1E-04*ABS(AREA12)
  3539.      *  .OR.ERRO12.LT.0.99E+00*ERRMAX) GO TO 10
  3540.         IF(EXTRAP) IROFF2 = IROFF2+1
  3541.         IF(.NOT.EXTRAP) IROFF1 = IROFF1+1
  3542.    10   IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF3 = IROFF3+1
  3543.    15   RLIST(MAXERR) = AREA1
  3544.         RLIST(LAST) = AREA2
  3545.         ERRBND = AMAX1(EPSABS,EPSREL*ABS(AREA))
  3546. C
  3547. C           TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
  3548. C
  3549.         IF(IROFF1+IROFF2.GE.10.OR.IROFF3.GE.20) IER = 2
  3550.         IF(IROFF2.GE.5) IERRO = 3
  3551. C
  3552. C           SET ERROR FLAG IN THE CASE THAT THE NUMBER OF SUBINTERVALS
  3553. C           EQUALS LIMIT.
  3554. C
  3555.         IF(LAST.EQ.LIMIT) IER = 1
  3556. C
  3557. C           SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
  3558. C           AT A POINT OF THE INTEGRATION RANGE.
  3559. C
  3560.         IF(AMAX1(ABS(A1),ABS(B2)).LE.(0.1E+01+0.1E+03*EPMACH)*
  3561.      *  (ABS(A2)+0.1E+04*UFLOW)) IER = 4
  3562. C
  3563. C           APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
  3564. C
  3565.         IF(ERROR2.GT.ERROR1) GO TO 20
  3566.         ALIST(LAST) = A2
  3567.         BLIST(MAXERR) = B1
  3568.         BLIST(LAST) = B2
  3569.         ELIST(MAXERR) = ERROR1
  3570.         ELIST(LAST) = ERROR2
  3571.         CALL QXCPY(VALP(1,MAXERR),VP1,LP1)
  3572.         LP(MAXERR)=LP1
  3573.         CALL QXCPY(VALN(1,MAXERR),VN1,LN1)
  3574.         LN(MAXERR)=LN1
  3575.         CALL QXCPY(VALP(1,LAST),VP2,LP2)
  3576.         LP(LAST)=LP2
  3577.         CALL QXCPY(VALN(1,LAST),VN2,LN2)
  3578.         LN(LAST)=LN2
  3579.         GO TO 30
  3580.    20   ALIST(MAXERR) = A2
  3581.         ALIST(LAST) = A1
  3582.         BLIST(LAST) = B1
  3583.         RLIST(MAXERR) = AREA2
  3584.         RLIST(LAST) = AREA1
  3585.         ELIST(MAXERR) = ERROR2
  3586.         ELIST(LAST) = ERROR1
  3587.         CALL QXCPY(VALP(1,MAXERR),VP2,LP2)
  3588.         LP(MAXERR)=LP2
  3589.         CALL QXCPY(VALN(1,MAXERR),VN2,LN2)
  3590.         LN(MAXERR)=LN2
  3591.         CALL QXCPY(VALP(1,LAST),VP1,LP1)
  3592.         LP(LAST)=LP1
  3593.         CALL QXCPY(VALN(1,LAST),VN1,LN1)
  3594.         LN(LAST)=LN1
  3595. C
  3596. C           CALL SUBROUTINE QPSRT TO MAINTAIN THE DESCENDING ORDERING
  3597. C           IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
  3598. C           WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
  3599. C
  3600.    30   CALL QPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX)
  3601. C ***JUMP OUT OF DO-LOOP
  3602.         IF(ERRSUM.LE.ERRBND) GO TO 115
  3603. C ***JUMP OUT OF DO-LOOP
  3604.         IF(IER.NE.0) GO TO 100
  3605.         IF(LAST.EQ.2) GO TO 80
  3606.         IF(NOEXT) GO TO 90
  3607.         ERLARG = ERLARG-ERLAST
  3608.         IF(ABS(B1-A1).GT.SMALL) ERLARG = ERLARG+ERRO12
  3609.         IF(EXTRAP) GO TO 40
  3610. C
  3611. C           TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE
  3612. C           SMALLEST INTERVAL.
  3613. C
  3614.         IF(ABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90
  3615.         EXTRAP = .TRUE.
  3616.         NRMAX = 2
  3617. C
  3618. C           THE BOUND 0.3*ERTEST HAS BEEN INTRODUCED TO PERFORM A
  3619. C           MORE CAUTIOUS EXTRAPOLATION THAN IN THE ORIGINAL DQAGSE
  3620. C           ROUTINE
  3621. C
  3622.    40   IF(IERRO.EQ.3.OR.ERLARG.LE.0.3E0*ERTEST) GO TO 60
  3623. C
  3624. C           THE SMALLEST INTERVAL HAS THE LARGEST ERROR.
  3625. C           BEFORE BISECTING DECREASE THE SUM OF THE ERRORS OVER THE
  3626. C           LARGER INTERVALS (ERLARG) AND PERFORM EXTRAPOLATION.
  3627. C
  3628.         ID = NRMAX
  3629.         JUPBND = LAST
  3630.         IF(LAST.GT.(2+LIMIT/2)) JUPBND = LIMIT+3-LAST
  3631.         DO 50 K = ID,JUPBND
  3632.           MAXERR = IORD(NRMAX)
  3633.           ERRMAX = ELIST(MAXERR)
  3634. C ***JUMP OUT OF DO-LOOP
  3635.           IF(ABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90
  3636.           NRMAX = NRMAX+1
  3637.    50   CONTINUE
  3638. C
  3639. C           PERFORM EXTRAPOLATION.
  3640. C
  3641.    60   NUMRL2 = NUMRL2+1
  3642.         RLIST2(NUMRL2) = AREA
  3643.         CALL QELG(NUMRL2,RLIST2,RESEPS,ABSEPS,RES3LA,NRES)
  3644.         KTMIN = KTMIN+1
  3645.         IF(KTMIN.GT.5.AND.ABSERR.LT.0.1E-02*ERRSUM) IER = 5
  3646.         IF(ABSEPS.GE.ABSERR) GO TO 70
  3647.         KTMIN = 0
  3648.         ABSERR = ABSEPS
  3649.         RESULT = RESEPS
  3650.         CORREC = ERLARG
  3651.         ERTEST = AMAX1(EPSABS,EPSREL*ABS(RESEPS))
  3652. C ***JUMP OUT OF DO-LOOP
  3653.         IF(ABSERR.LE.ERTEST) GO TO 100
  3654. C
  3655. C           PREPARE BISECTION OF THE SMALLEST INTERVAL.
  3656. C
  3657.    70   IF(NUMRL2.EQ.1) NOEXT = .TRUE.
  3658.         IF(IER.EQ.5) GO TO 100
  3659.         MAXERR = IORD(1)
  3660.         ERRMAX = ELIST(MAXERR)
  3661.         NRMAX = 1
  3662.         EXTRAP = .FALSE.
  3663.         SMALL = SMALL*0.5E+00
  3664.         ERLARG = ERRSUM
  3665.         GO TO 90
  3666.    80   SMALL = ABS(B-A)*0.375E+00
  3667.         ERLARG = ERRSUM
  3668.         ERTEST = ERRBND
  3669.         RLIST2(2) = AREA
  3670.    90 CONTINUE
  3671. C
  3672. C           SET FINAL RESULT AND ERROR ESTIMATE.
  3673. C           ------------------------------------
  3674. C
  3675.   100 IF(ABSERR.EQ.OFLOW) GO TO 115
  3676.       IF(IER+IERRO.EQ.0) GO TO 110
  3677.       IF(IERRO.EQ.3) ABSERR = ABSERR+CORREC
  3678.       IF(IER.EQ.0) IER = 3
  3679.       IF(RESULT.NE.0.0E+00.AND.AREA.NE.0.0E+00) GO TO 105
  3680.       IF(ABSERR.GT.ERRSUM) GO TO 115
  3681.       IF(AREA.EQ.0.0E+00) GO TO 130
  3682.       GO TO 110
  3683.   105 IF(ABSERR/ABS(RESULT).GT.ERRSUM/ABS(AREA)) GO TO 115
  3684. C
  3685. C           TEST ON DIVERGENCE.
  3686. C
  3687.   110 IF(KSGN.EQ.(-1).AND.AMAX1(ABS(RESULT),ABS(AREA)).LE.
  3688.      * DEFABS*0.1E-01) GO TO 130
  3689.       IF(0.1E-01.GT.(RESULT/AREA).OR.(RESULT/AREA).GT.0.1E+03
  3690.      * .OR.ERRSUM.GT.ABS(AREA)) IER = 6
  3691.       GO TO 130
  3692. C
  3693. C           COMPUTE GLOBAL INTEGRAL SUM.
  3694. C
  3695.   115 RESULT = 0.0E+00
  3696.       DO 120 K = 1,LAST
  3697.          RESULT = RESULT+RLIST(K)
  3698.   120 CONTINUE
  3699.       ABSERR = ERRSUM
  3700.   130 IF(IER.GT.2) IER = IER-1
  3701.   999 RETURN
  3702.       END
  3703. *DECK QXG
  3704.       SUBROUTINE QXG(F,A,B,EPSABS,EPSREL,KEY,RESULT,ABSERR,IER,
  3705.      1   LIMIT,LENIW,LENW,LAST,IWORK,WORK)
  3706. C
  3707. C            THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
  3708. C            DEFINITE INTEGRAL  I = INTEGRAL OF F OVER (A,B),
  3709. C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
  3710. C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
  3711. C
  3712. C        PARAMETERS
  3713. C         ON ENTRY
  3714. C            F      - REAL
  3715. C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  3716. C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  3717. C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  3718. C
  3719. C            A      - REAL
  3720. C                     LOWER LIMIT OF INTEGRATION
  3721. C
  3722. C            B      - REAL
  3723. C                     UPPER LIMIT OF INTEGRATION
  3724. C
  3725. C            EPSABS - REAL
  3726. C                     ABSOLUTE ACCURACY REQUESTED
  3727. C            EPSREL - REAL
  3728. C                     RELATIVE ACCURACY REQUESTED
  3729. C                     IF  EPSABS.LE.0
  3730. C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28),
  3731. C                     THE ROUTINE WILL END WITH IER = 6.
  3732. C
  3733. C            KEY    - INTEGER
  3734. C                     KEY FOR CHOICE OF LOCAL INTEGRATION RULE
  3735. C                     RMS FORMULAS ARE USED WITH
  3736. C                      13 - 19               POINTS IF KEY.LT.1,
  3737. C                      13 - 19 - (27)        POINTS IF KEY = 1,
  3738. C                      13 - 19 - (27) - (41) POINTS IF KEY = 2,
  3739. C                           19 -  27  - (41) POINTS IF KEY = 3,
  3740. C                                 27  -  41  POINTS IF KEY.GT.3.
  3741. C
  3742. C                         (RULES) USED IF THE FUNCTION APPEARS
  3743. C                         ENOUGH REGULAR IN THE SUBINTERVAL
  3744. C
  3745. C         ON RETURN
  3746. C            RESULT - REAL
  3747. C                     APPROXIMATION TO THE INTEGRAL
  3748. C
  3749. C            ABSERR - REAL
  3750. C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
  3751. C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
  3752. C
  3753. C            IER    - INTEGER
  3754. C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
  3755. C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
  3756. C                             ACCURACY HAS BEEN ACHIEVED.
  3757. C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE
  3758. C                             THE ESTIMATES FOR RESULT AND ERROR ARE
  3759. C                             LESS RELIABLE. IT IS ASSUMED THAT THE
  3760. C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
  3761. C            ERROR MESSAGES
  3762. C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
  3763. C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB-
  3764. C                             DIVISIONS BY INCREASING THE VALUE OF LIMIT
  3765. C                             (AND TAKING THE ACCORDING DIMENSION
  3766. C                             ADJUSTMENTS INTO ACCOUNT. HOWEVER, IF
  3767. C                             THIS YIELDS NO IMPROVEMENT IT IS ADVISED
  3768. C                             TO ANALYZE THE INTEGRAND IN ORDER TO
  3769. C                             DETERMINE THE INTEGRATION DIFFICULTIES. IF
  3770. C                             THE POSITION OF A LOCAL DIFFICULTY CAN BE
  3771. C                             DETERMINED (E.G. SINGULARITY,
  3772. C                             DISCONTINUITY WITHIN THE INTERVAL) ONE
  3773. C                             WILL PROBABLY GAIN FROM SPLITTING UP THE
  3774. C                             INTERVAL AT THIS POINT AND CALLING THE
  3775. C                             INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
  3776. C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
  3777. C                             SHOULD BE USED, WHICH IS DESIGNED FOR
  3778. C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.
  3779. C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS DETEC-
  3780. C                             TED, WHICH PREVENTS THE REQUESTED
  3781. C                             TOLERANCE FROM BEING ACHIEVED.
  3782. C                             THE ERROR MAY BE UNDER-ESTIMATED.
  3783. C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR
  3784. C                             OCCURS AT SOME POINTS OF THE INTEGRATION
  3785. C                             INTERVAL.
  3786. C                         = 4 THE ALGORITHM DOES NOT CONVERGE.
  3787. C                             ROUNDOFF ERROR IS DETECTED IN THE
  3788. C                             EXTRAPOLATION TABLE. IT IS PRESUMED THAT
  3789. C                             THE REQUESTED TOLERANCE CANNOT BE
  3790. C                             ACHIEVED, AND THAT THE RETURNED RESULT IS
  3791. C                             THE BEST WHICH CAN BE OBTAINED.
  3792. C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
  3793. C                             SLOWLY CONVERGENT. IT MUST BE NOTED THAT
  3794. C                             DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
  3795. C                             OF IER.
  3796. C                         = 6 THE INPUT IS INVALID, BECAUSE
  3797. C                             (EPSABS.LE.0 AND
  3798. C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28)
  3799. C                             OR LIMIT.LT.1 OR LENW.LT.LIMIT*46 OR
  3800. C                                LENIW .LT. LIMIT*3
  3801. C                             RESULT, ABSERR, LAST ARE SET TO
  3802. C                             ZERO.EXCEPT WHEN LIMIT OR LENW OR LENIW
  3803. C                             IS INVALID, IWORK(1), WORK(LIMIT*2+1) AND
  3804. C                             WORK(LIMIT*3+1) ARE SET TO ZERO, WORK(1)
  3805. C                             IS SET TO A AND WORK(LIMIT+1) TO B.
  3806. C
  3807. C         DIMENSIONING PARAMETERS
  3808. C            LIMIT - INTEGER
  3809. C                    LIMIT DETERMINES THE MAXIMUM NUMBER OF SUBINTERVALS
  3810. C                    IN THE PARTITION OF THE GIVEN INTEGRATION INTERVAL
  3811. C                    (A,B), LIMIT.GE.1.
  3812. C                    IF LIMIT.LT.1, THE ROUTINE WILL END WITH IER = 6.
  3813. C
  3814. C            LENW  - INTEGER
  3815. C                    DIMENSIONING PARAMETER FOR WORK
  3816. C                    LENW MUST BE AT LEAST LIMIT*46
  3817. C                    IF LENW.LT.LIMIT*46, THE ROUTINE WILL END
  3818. C                    WITH IER = 6.
  3819. C
  3820. C            LENIW - INTEGER
  3821. C                    DIMENSIONING PARAMETER FOR IWORK
  3822. C                    LENIW MUST BE AT LEAST LIMIT*3
  3823. C                    IF LENW.LT.LIMIT*3, THE ROUTINE WILL END
  3824. C                    WITH IER = 6.
  3825. C
  3826. C            LAST  - INTEGER
  3827. C                    ON RETURN, LAST EQUALS THE NUMBER OF SUBINTERVALS
  3828. C                    PRODUCED IN THE SUBDIVISION PROCESS, DETERMINES THE
  3829. C                    NUMBER OF SIGNIFICANT ELEMENTS ACTUALLY IN THE WORK
  3830. C                    ARRAYS.
  3831. C
  3832. C         WORK ARRAYS
  3833. C            IWORK - INTEGER
  3834. C                    VECTOR OF DIMENSION AT LEAST 3*LIMIT, THE FIRST K
  3835. C                    ELEMENTS OF WHICH CONTAIN POINTERS
  3836. C                    TO THE ERROR ESTIMATES OVER THE SUBINTERVALS
  3837. C                    SUCH THAT WORK(LIMIT*3+IWORK(1)),... ,
  3838. C                    WORK(LIMIT*3+IWORK(K)) FORM A DECREASING
  3839. C                    SEQUENCE, WITH K = LAST IF LAST.LE.(LIMIT/2+2),
  3840. C                    AND K = LIMIT+1-LAST OTHERWISE
  3841. C
  3842. C            WORK  - REAL
  3843. C                    VECTOR OF DIMENSION AT LEAST LENW
  3844. C                    ON RETURN
  3845. C                    WORK(1), ..., WORK(LAST) CONTAIN THE LEFT
  3846. C                     END-POINTS OF THE SUBINTERVALS IN THE
  3847. C                     PARTITION OF (A,B),
  3848. C                    WORK(LIMIT+1), ..., WORK(LIMIT+LAST) CONTAIN
  3849. C                     THE RIGHT END-POINTS,
  3850. C                    WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) CONTAIN
  3851. C                     THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS,
  3852. C                    WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST)
  3853. C                     CONTAIN THE ERROR ESTIMATES.
  3854. C                    WORK(LIMIT*4+1), ... IS THE AREA RESERVED TO STORE
  3855. C                     FUNCTIONAL VALUES.
  3856. C***ROUTINES CALLED  QXGE,XERROR
  3857. C
  3858. C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
  3859. C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
  3860. C
  3861. C
  3862.       REAL A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK
  3863.       INTEGER KEY,IER,IWORK,LAST,LENIW,LENW,LIMIT,LVL,L1,L2,L3,L4,L5
  3864. C
  3865.       DIMENSION IWORK(LENIW),WORK(LENW)
  3866. C
  3867.       EXTERNAL F
  3868. C
  3869. C         CHECK VALIDITY OF LIMIT,LENIW AND LENW.
  3870. C
  3871. C***FIRST EXECUTABLE STATEMENT  QXG
  3872.       IER = 6
  3873.       LAST = 0
  3874.       RESULT = 0.0E+00
  3875.       ABSERR = 0.0E+00
  3876.       IF(LIMIT.LT.1.OR.LENIW.LT.LIMIT*3.OR.LENW.LT.LIMIT*46) GO TO 10
  3877. C
  3878. C         PREPARE CALL FOR QXGE.
  3879. C
  3880.       L1 = LIMIT+1
  3881.       L2 = LIMIT+L1
  3882.       L3 = LIMIT+L2
  3883.       L4 = LIMIT+L3
  3884.       L5 = 21*LIMIT+L3
  3885. C
  3886.       CALL QXGE(F,A,B,EPSABS,EPSREL,KEY,LIMIT,RESULT,ABSERR,
  3887.      *  IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST,
  3888.      *  WORK(L4),WORK(L5),IWORK(L1),IWORK(L2))
  3889. C
  3890. C         CALL ERROR HANDLER IF NECESSARY.
  3891. C
  3892.       LVL = 0
  3893. 10    IF(IER.EQ.6) LVL = 1
  3894. C
  3895. C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
  3896. C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
  3897. C
  3898. C     IF(IER.NE.0)CALL XERROR('ABNORMAL RETURN FROM QXG   ',28,IER,LVL)
  3899.       RETURN
  3900.       END
  3901. *DECK QXGE
  3902.       SUBROUTINE QXGE(F,A,B,EPSABS,EPSREL,KEY,LIMIT,RESULT,ABSERR,
  3903.      *   IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST,VALP,VALN,LP,LN)
  3904. C
  3905. C
  3906. C            THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
  3907. C            DEFINITE INTEGRAL   I = INTEGRAL OF F OVER (A,B),
  3908. C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
  3909. C            ABS(I-RESLT).LE.MAX(EPSABS,EPSREL*ABS(I)).
  3910. C
  3911. C        PARAMETERS
  3912. C         ON ENTRY
  3913. C            F      - REAL
  3914. C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  3915. C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  3916. C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  3917. C
  3918. C            A      - REAL
  3919. C                     LOWER LIMIT OF INTEGRATION
  3920. C
  3921. C            B      - REAL
  3922. C                     UPPER LIMIT OF INTEGRATION
  3923. C
  3924. C            EPSABS - REAL
  3925. C                     ABSOLUTE ACCURACY REQUESTED
  3926. C
  3927. C            EPSREL - REAL
  3928. C                     RELATIVE ACCURACY REQUESTED
  3929. C                     IF  EPSABS.LE.0
  3930. C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28),
  3931. C                     THE ROUTINE WILL END WITH IER = 6.
  3932. C
  3933. C            KEY    - INTEGER
  3934. C                     KEY FOR CHOICE OF LOCAL INTEGRATION RULE
  3935. C                     RMS FORMULAS ARE USED WITH
  3936. C                      13 - 19               POINTS IF KEY.LT.1,
  3937. C                      13 - 19 - (27)        POINTS IF KEY = 1,
  3938. C                      13 - 19 - (27) - (41) POINTS IF KEY = 2,
  3939. C                           19 -  27  - (41) POINTS IF KEY = 3,
  3940. C                                 27  -  41  POINTS IF KEY.GT.3.
  3941. C
  3942. C                         (RULES) USED IF THE FUNCTION APPEARS
  3943. C                         ENOUGH REGULAR IN THE SUBINTERVAL
  3944. C
  3945. C            LIMIT  - INTEGER
  3946. C                     GIVES AN UPPERBOUND ON THE NUMBER OF SUBINTERVALS
  3947. C                     IN THE PARTITION OF (A,B), LIMIT.GE.1.
  3948. C
  3949. C         ON RETURN
  3950. C            RESULT - REAL
  3951. C                     APPROXIMATION TO THE INTEGRAL
  3952. C
  3953. C            ABSERR - REAL
  3954. C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
  3955. C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
  3956. C
  3957. C            IER    - INTEGER
  3958. C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
  3959. C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
  3960. C                             ACCURACY HAS BEEN ACHIEVED.
  3961. C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE
  3962. C                             THE ESTIMATES FOR RESULT AND ERROR ARE
  3963. C                             LESS RELIABLE. IT IS ASSUMED THAT THE
  3964. C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
  3965. C            ERROR MESSAGES
  3966. C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
  3967. C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE
  3968. C                             SUBDIVISIONS BY INCREASING THE VALUE
  3969. C                             OF LIMIT.
  3970. C                             HOWEVER, IF THIS YIELDS NO IMPROVEMENT IT
  3971. C                             IS RATHER ADVISED TO ANALYZE THE INTEGRAND
  3972. C                             IN ORDER TO DETERMINE THE INTEGRATION
  3973. C                             DIFFICULTIES. IF THE POSITION OF A LOCAL
  3974. C                             DIFFICULTY CAN BE DETERMINED(E.G.
  3975. C                             SINGULARITY, DISCONTINUITY WITHIN THE
  3976. C                             INTERVAL) ONE WILL PROBABLY GAIN FROM
  3977. C                             SPLITTING UP THE INTERVAL AT THIS POINT
  3978. C                             AND CALLING THE INTEGRATOR ON THE
  3979. C                             SUBRANGES. IF POSSIBLE, AN APPROPRIATE
  3980. C                             SPECIAL-PURPOSE INTEGRATOR SHOULD BE USED
  3981. C                             WHICH IS DESIGNED FOR HANDLING THE TYPE OF
  3982. C                             DIFFICULTY INVOLVED.
  3983. C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS
  3984. C                             DETECTED, WHICH PREVENTS THE REQUESTED
  3985. C                             TOLERANCE FROM BEING ACHIEVED.
  3986. C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS
  3987. C                             AT SOME POINTS OF THE INTEGRATION
  3988. C                             INTERVAL.
  3989. C                         = 6 THE INPUT IS INVALID, BECAUSE
  3990. C                             (EPSABS.LE.0 AND
  3991. C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28),
  3992. C                             RESULT, ABSERR, LAST, RLIST(1) ,
  3993. C                             ELIST(1) AND IORD(1) ARE SET TO ZERO.
  3994. C                             ALIST(1) AND BLIST(1) ARE SET TO A AND B
  3995. C                             RESPECTIVELY.
  3996. C
  3997. C            ALIST   - REAL
  3998. C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  3999. C                       LAST  ELEMENTS OF WHICH ARE THE LEFT
  4000. C                      END POINTS OF THE SUBINTERVALS IN THE PARTITION
  4001. C                      OF THE GIVEN INTEGRATION RANGE (A,B)
  4002. C
  4003. C            BLIST   - REAL
  4004. C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  4005. C                       LAST  ELEMENTS OF WHICH ARE THE RIGHT
  4006. C                      END POINTS OF THE SUBINTERVALS IN THE PARTITION
  4007. C                      OF THE GIVEN INTEGRATION RANGE (A,B)
  4008. C
  4009. C            RLIST   - REAL
  4010. C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  4011. C                       LAST  ELEMENTS OF WHICH ARE THE
  4012. C                      INTEGRAL APPROXIMATIONS ON THE SUBINTERVALS
  4013. C
  4014. C            ELIST   - REAL
  4015. C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  4016. C                       LAST  ELEMENTS OF WHICH ARE THE MODULI OF THE
  4017. C                      ABSOLUTE ERROR ESTIMATES ON THE SUBINTERVALS
  4018. C
  4019. C            IORD    - INTEGER
  4020. C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST K
  4021. C                      ELEMENTS OF WHICH ARE POINTERS TO THE
  4022. C                      ERROR ESTIMATES OVER THE SUBINTERVALS,
  4023. C                      SUCH THAT ELIST(IORD(1)), ...,
  4024. C                      ELIST(IORD(K)) FORM A DECREASING SEQUENCE,
  4025. C                      WITH K = LAST IF LAST.LE.(LIMIT/2+2), AND
  4026. C                      K = LIMIT+1-LAST OTHERWISE
  4027. C
  4028. C            LAST    - INTEGER
  4029. C                      NUMBER OF SUBINTERVALS ACTUALLY PRODUCED IN THE
  4030. C                      SUBDIVISION PROCESS
  4031. C
  4032. C            VALP   - REAL
  4033. C            VALN     ARRAYS OF DIMENSION AT LEAST (21,LIMIT) USED TO
  4034. C                     SAVE THE FUNCTIONAL VALUES
  4035. C
  4036. C            LP     - INTEGER
  4037. C            LN       VECTORS OF DIMENSION AT LEAST LIMIT, USED TO
  4038. C                     STORE THE ACTUAL NUMBER OF FUNCTIONAL VALUES
  4039. C                     SAVED IN THE CORRESPONDING COLUMN
  4040. C                     OF VALP,VALN
  4041. C
  4042. C***ROUTINES CALLED  F,R1MACH,QXLQM,QXRRD,QPSRT,QXCPY
  4043. C
  4044.       REAL A,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,A2,B,
  4045.      *  BLIST,B1,B2,ABS,DEFABS,DEFAB1,DEFAB2,AMAX1,R1MACH,ELIST,EPMACH,
  4046.      *  EPSABS,EPSREL,ERRBND,ERRMAX,ERROR1,ERROR2,ERRO12,ERRSUM,F,
  4047.      *  RESABS,RESULT,RLIST,UFLOW,VALP,VALN,VP1,VP2,VN1,VN2
  4048.       INTEGER KEY,IER,IORD,IROFF1,IROFF2,K,LAST,LIMIT,MAXERR,
  4049.      *  NRMAX,LP,LN,LP1,LP2,LN1,LN2
  4050. C
  4051.       DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT),
  4052.      *  RLIST(LIMIT),VALP(21,LIMIT),VALN(21,LIMIT),LP(LIMIT),LN(LIMIT),
  4053.      * VP1(21),VP2(21),VN1(21),VN2(21)
  4054. C
  4055.       EXTERNAL F
  4056. C
  4057. C            MACHINE DEPENDENT CONSTANTS
  4058. C            ---------------------------
  4059. C
  4060. C          EPMACH IS THE LARGEST RELATIVE SPACING.
  4061. C          UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
  4062. C
  4063. C***FIRST EXECUTABLE STATEMENT  QXGE
  4064.       EPMACH = R1MACH(4)
  4065.       UFLOW = R1MACH(1)
  4066. C
  4067. C           TEST ON VALIDITY OF PARAMETERS
  4068. C           ------------------------------
  4069. C
  4070.       IER = 0
  4071.       LAST = 0
  4072.       RESULT = 0.0E+00
  4073.       ABSERR = 0.0E+00
  4074.       ALIST(1) = A
  4075.       BLIST(1) = B
  4076.       RLIST(1) = 0.0E+00
  4077.       ELIST(1) = 0.0E+00
  4078.       IORD(1) = 0
  4079.       IF(EPSABS.LE.0.0E+00.AND.
  4080.      *  EPSREL.LT.AMAX1(0.5E+02*EPMACH,0.5E-28)) IER = 6
  4081.       IF(IER.EQ.6) GO TO 999
  4082. C
  4083. C           FIRST APPROXIMATION TO THE INTEGRAL
  4084. C           -----------------------------------
  4085. C
  4086.       LP(1)=1
  4087.       LN(1)=1
  4088.       VALP(1,1)=F((A+B)*0.5E0)
  4089.       VALN(1,1)=VALP(1,1)
  4090.       CALL QXLQM(F,A,B,RESULT,ABSERR,DEFABS,RESABS,
  4091.      *           VALP(1,1),VALN(1,1),LP(1),LN(1),KEY)
  4092.       LAST = 1
  4093.       RLIST(1) = RESULT
  4094.       ELIST(1) = ABSERR
  4095.       IORD(1) = 1
  4096. C
  4097. C           TEST ON ACCURACY.
  4098. C
  4099.       ERRBND = AMAX1(EPSABS,EPSREL*ABS(RESULT))
  4100.       IF(ABSERR.LE.0.5E+02*EPMACH*DEFABS.AND.ABSERR.GT.ERRBND) IER = 2
  4101.       IF(LIMIT.EQ.1) IER = 1
  4102.       IF(IER.NE.0.OR.(ABSERR.LE.ERRBND.AND.ABSERR.NE.RESABS)
  4103.      *  .OR.ABSERR.EQ.0.0E+00) GO TO 999
  4104. C
  4105. C           INITIALIZATION
  4106. C           --------------
  4107. C
  4108. C
  4109.       ERRMAX = ABSERR
  4110.       MAXERR = 1
  4111.       AREA = RESULT
  4112.       ERRSUM = ABSERR
  4113.       NRMAX = 1
  4114.       IROFF1 = 0
  4115.       IROFF2 = 0
  4116. C
  4117. C           MAIN DO-LOOP
  4118. C           ------------
  4119. C
  4120.       DO 30 LAST = 2,LIMIT
  4121. C
  4122. C           BISECT THE SUBINTERVAL WITH THE LARGEST ERROR ESTIMATE.
  4123. C
  4124.         A1 = ALIST(MAXERR)
  4125.         B1 = 0.5E+00*(ALIST(MAXERR)+BLIST(MAXERR))
  4126.         A2 = B1
  4127.         B2 = BLIST(MAXERR)
  4128.         CALL QXRRD(F,VALN(1,MAXERR),LN(MAXERR),B1,A1,VN1,VP1,LN1,LP1)
  4129.         CALL QXLQM(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1,VP1,VN1,LP1,LN1,
  4130.      *              KEY)
  4131.         CALL QXRRD(F,VALP(1,MAXERR),LP(MAXERR),A2,B2,VP2,VN2,LP2,LN2)
  4132.         CALL QXLQM(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2,VP2,VN2,LP2,LN2,
  4133.      *              KEY)
  4134. C
  4135. C           IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
  4136. C           AND ERROR AND TEST FOR ACCURACY.
  4137. C
  4138.         AREA12 = AREA1+AREA2
  4139.         ERRO12 = ERROR1+ERROR2
  4140.         ERRSUM = ERRSUM+ERRO12-ERRMAX
  4141.         AREA = AREA+AREA12-RLIST(MAXERR)
  4142.         IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2) GO TO 5
  4143.         IF(ABS(RLIST(MAXERR)-AREA12).LE.0.1E-04*ABS(AREA12)
  4144.      *  .AND.ERRO12.GE.0.99E+00*ERRMAX) IROFF1 = IROFF1+1
  4145.         IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF2 = IROFF2+1
  4146.     5   RLIST(MAXERR) = AREA1
  4147.         RLIST(LAST) = AREA2
  4148.         ERRBND = AMAX1(EPSABS,EPSREL*ABS(AREA))
  4149.         IF(ERRSUM.LE.ERRBND) GO TO 8
  4150. C
  4151. C           TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
  4152. C
  4153.         IF(IROFF1.GE.6.OR.IROFF2.GE.20) IER = 2
  4154. C
  4155. C           SET ERROR FLAG IN THE CASE THAT THE NUMBER OF SUBINTERVALS
  4156. C           EQUALS LIMIT.
  4157. C
  4158.         IF(LAST.EQ.LIMIT) IER = 1
  4159. C
  4160. C           SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
  4161. C           AT A POINT OF THE INTEGRATION RANGE.
  4162. C
  4163.         IF(AMAX1(ABS(A1),ABS(B2)).LE.(0.1E+01+0.1E+03*
  4164.      *  EPMACH)*(ABS(A2)+0.1E+04*UFLOW)) IER = 3
  4165. C
  4166. C           APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
  4167. C
  4168.     8   IF(ERROR2.GT.ERROR1) GO TO 10
  4169.         ALIST(LAST) = A2
  4170.         BLIST(MAXERR) = B1
  4171.         BLIST(LAST) = B2
  4172.         ELIST(MAXERR) = ERROR1
  4173.         ELIST(LAST) = ERROR2
  4174.         CALL QXCPY(VALP(1,MAXERR),VP1,LP1)
  4175.         LP(MAXERR)=LP1
  4176.         CALL QXCPY(VALN(1,MAXERR),VN1,LN1)
  4177.         LN(MAXERR)=LN1
  4178.         CALL QXCPY(VALP(1,LAST),VP2,LP2)
  4179.         LP(LAST)=LP2
  4180.         CALL QXCPY(VALN(1,LAST),VN2,LN2)
  4181.         LN(LAST)=LN2
  4182.         GO TO 20
  4183.    10   ALIST(MAXERR) = A2
  4184.         ALIST(LAST) = A1
  4185.         BLIST(LAST) = B1
  4186.         RLIST(MAXERR) = AREA2
  4187.         RLIST(LAST) = AREA1
  4188.         ELIST(MAXERR) = ERROR2
  4189.         ELIST(LAST) = ERROR1
  4190.         CALL QXCPY(VALP(1,MAXERR),VP2,LP2)
  4191.         LP(MAXERR)=LP2
  4192.         CALL QXCPY(VALN(1,MAXERR),VN2,LN2)
  4193.         LN(MAXERR)=LN2
  4194.         CALL QXCPY(VALP(1,LAST),VP1,LP1)
  4195.         LP(LAST)=LP1
  4196.         CALL QXCPY(VALN(1,LAST),VN1,LN1)
  4197.         LN(LAST)=LN1
  4198. C
  4199. C           CALL SUBROUTINE QPSRT TO MAINTAIN THE DESCENDING ORDERING
  4200. C           IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
  4201. C           WITH THE LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
  4202. C
  4203.    20   CALL QPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX)
  4204. C ***JUMP OUT OF DO-LOOP
  4205.         IF(IER.NE.0.OR.ERRSUM.LE.ERRBND) GO TO 40
  4206.    30 CONTINUE
  4207. C
  4208. C           COMPUTE FINAL RESULT.
  4209. C           ---------------------
  4210. C
  4211.    40 RESULT = 0.0E+00
  4212.       DO 50 K=1,LAST
  4213.         RESULT = RESULT+RLIST(K)
  4214.    50 CONTINUE
  4215.       ABSERR = ERRSUM
  4216.   999 RETURN
  4217.       END
  4218. *DECK QXLQM
  4219.       SUBROUTINE QXLQM(F,A,B,RESULT,ABSERR,RESABS,RESASC,VR,VS,LR,LS,
  4220.      *                  KEY)
  4221. C
  4222. C            TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR
  4223. C                           ESTIMATE
  4224. C                       J = INTEGRAL OF ABS(F) OVER (A,B)
  4225. C
  4226. C           PARAMETERS
  4227. C            ON ENTRY
  4228. C              F      - REAL
  4229. C                       FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  4230. C                       FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  4231. C                       DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  4232. C
  4233. C              A      - REAL
  4234. C                       LOWER LIMIT OF INTEGRATION
  4235. C
  4236. C              B      - REAL
  4237. C                       UPPER LIMIT OF INTEGRATION
  4238. C
  4239. C              VR     - REAL
  4240. C                       VECTOR OF LENGTH LR CONTAINING THE
  4241. C                       SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS
  4242. C
  4243. C              VS     - REAL
  4244. C                       VECTOR OF LENGTH LS CONTAINING THE
  4245. C                       SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS
  4246. C
  4247. C              LR     - INTEGER
  4248. C              LS       NUMBER OF ELEMENTS IN
  4249. C                       VR,VS RESPECTIVELY
  4250. C
  4251. C            KEY    - INTEGER
  4252. C                     KEY FOR CHOICE OF LOCAL INTEGRATION RULE
  4253. C                     RMS FORMULAS ARE USED WITH
  4254. C                      13 - 19               POINTS IF KEY.LT.1,
  4255. C                      13 - 19 - (27)        POINTS IF KEY = 1,
  4256. C                      13 - 19 - (27) - (41) POINTS IF KEY = 2,
  4257. C                           19 -  27  - (41) POINTS IF KEY = 3,
  4258. C                                 27  -  41  POINTS IF KEY.GT.3.
  4259. C
  4260. C                         (RULES) USED IF THE FUNCTION APPEARS
  4261. C                         ENOUGH REGULAR
  4262. C
  4263. C            ON RETURN
  4264. C              RESULT - REAL
  4265. C                       APPROXIMATION TO THE INTEGRAL I
  4266. C
  4267. C              ABSERR - REAL
  4268. C                       ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
  4269. C                       WHICH SHOULD NOT EXCEED ABS(I-RESULT)
  4270. C
  4271. C              RESABS - REAL
  4272. C                       APPROXIMATION TO THE INTEGRAL J
  4273. C
  4274. C              RESASC - REAL
  4275. C                       APPROXIMATION TO THE INTEGRAL OF ABS(F-I/(B-A))
  4276. C                       OVER (A,B)
  4277. C
  4278. C              VR     - REAL
  4279. C                       VECTOR OF LENGTH LR CONTAINING THE
  4280. C                       SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS
  4281. C
  4282. C              VS     - REAL
  4283. C                       VECTOR OF LENGTH LS CONTAINING THE
  4284. C                       SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS
  4285. C
  4286. C              LR     - INTEGER
  4287. C              LS       NUMBER OF ELEMENTS IN
  4288. C                       VR,VS RESPECTIVELY
  4289. C
  4290. C***ROUTINES CALLED  R1MACH,QXRUL
  4291. C
  4292.       REAL F,A,B,RESULT,ABSERR,RESABS,RESASC,
  4293.      *          R1MACH,EPMACH,RESG,RESK,UFLOW,ERROLD,VR(21),VS(21)
  4294.       INTEGER K,K0,K1,K2,KEY,KEY1,LR,LS
  4295.       EXTERNAL F
  4296. C
  4297. C            MACHINE DEPENDENT CONSTANTS
  4298. C            ---------------------------
  4299. C
  4300. C          EPMACH IS THE LARGEST RELATIVE SPACING.
  4301. C          UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
  4302. C          ERROLD IS THE LARGEST MAGNITUDE.
  4303. C
  4304. C***FIRST EXECUTABLE STATEMENT QXLQM
  4305.       EPMACH = R1MACH(4)
  4306.       UFLOW  = R1MACH(1)
  4307.       ERROLD = R1MACH(2)
  4308. C
  4309.       KEY1 = MAX(KEY ,  0)
  4310.       KEY1 = MIN(KEY1,  4)
  4311.       K0   = MAX(KEY1-2,0)
  4312.       K1   = K0+1
  4313.       K2   = MIN(KEY1+1,3)
  4314. C
  4315.       CALL QXRUL(F,A,B,RESG,RESABS,RESASC,K0,K1,VR,VS,LR,LS)
  4316.       DO 99 K=K1,K2
  4317.         CALL QXRUL(F,A,B,RESK,RESABS,RESASC,K,K1,VR,VS,LR,LS)
  4318.         RESULT=RESK
  4319.         ABSERR = ABS((RESK-RESG))
  4320.         IF(RESASC.NE.0.0E+00.AND.ABSERR.NE.0.0E+00)
  4321.      *  ABSERR = RESASC*AMIN1(1.E0,(200E0*ABSERR/RESASC)**1.5E+00)
  4322.         IF(RESABS.GT.UFLOW/(10E0*EPMACH)) ABSERR = AMAX1
  4323.      *  ((EPMACH*10E0)*RESABS,ABSERR)
  4324.         RESG=RESK
  4325.         IF(ABSERR.GT.ERROLD*0.16)GOTO 3000
  4326.         IF(ABSERR.LT.1000*EPMACH*RESABS)GOTO 3000
  4327.         ERROLD=ABSERR
  4328. 99    CONTINUE
  4329. 3000  CONTINUE
  4330.       RETURN
  4331.       END
  4332. *DECK QXRUL
  4333.       SUBROUTINE QXRUL(F,XL,XU,Y,YA,YM,KE,K1,FV1,FV2,L1,L2)
  4334. C
  4335. C            TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR
  4336. C                           ESTIMATE
  4337. C            AND CONDITIONALLY COMPUTE
  4338. C                       J = INTEGRAL OF ABS(F) OVER (A,B)
  4339. C                       BY USING AN  RMS RULE
  4340. C           PARAMETERS
  4341. C            ON ENTRY
  4342. C              F      - REAL
  4343. C                       FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  4344. C                       FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  4345. C                       DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  4346. C
  4347. C              XL     - REAL
  4348. C                       LOWER LIMIT OF INTEGRATION
  4349. C
  4350. C              XU     - REAL
  4351. C                       UPPER LIMIT OF INTEGRATION
  4352. C
  4353. C
  4354. C              KE     - INTEGER
  4355. C                     KEY FOR CHOICE OF LOCAL INTEGRATION RULE
  4356. C                     AN RMS RULE IS USED WITH
  4357. C                         13      POINTS IF KE  = 2,
  4358. C                         19      POINTS IF KE  = 3,
  4359. C                         27      POINTS IF KE  = 4,
  4360. C                         42      POINTS IF KE  = 5
  4361. C
  4362. C              K1     INTEGER
  4363. C                     VALUE OF KEY FOR WHICH THE ADDITIONAL ESTIMATES
  4364. C                     YA, YM ARE TO BE COMPUTED
  4365. C
  4366. C              FV1    - REAL
  4367. C                       VECTOR CONTAINING L1
  4368. C                       SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS
  4369. C
  4370. C              FV2    - REAL
  4371. C                       VECTOR CONTAINING L2
  4372. C                       SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS
  4373. C
  4374. C              L1     - INTEGER
  4375. C              L2       NUMBER OF ELEMENTS IN FV1,FV2  RESPECTIVELY
  4376. C
  4377. C            ON RETURN
  4378. C              Y      - REAL
  4379. C                       APPROXIMATION TO THE INTEGRAL I
  4380. C                       RESULT IS COMPUTED BY APPLYING THE
  4381. C                       REQUESTED RMS RULE
  4382. C
  4383. C              YA     - REAL
  4384. C                       IF KEY = K1  APPROXIMATION TO THE INTEGRAL J
  4385. C                       ELSE UNCHANGED
  4386. C
  4387. C              YM     - REAL
  4388. C                       IF KEY = K1  APPROXIMATION TO THE INTEGRAL OF
  4389. C                                      ABS(F-I/(XU-XL)   OVER (XL,XU)
  4390. C                       ELSE UNCHANGED
  4391. C
  4392. C              FV1    - REAL
  4393. C                       VECTOR L1 CONTAINING L1
  4394. C                       SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS
  4395. C
  4396. C              FV2    - REAL
  4397. C                       VECTOR CONTAINING L2
  4398. C                       SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS
  4399. C
  4400. C              L1     - INTEGER
  4401. C              L2       NUMBER OF ELEMENTS IN FV1,FV2  RESPECTIVELY
  4402. C
  4403. C***ROUTINES CALLED  F
  4404. C
  4405.       REAL F,XL,XU,LDL,Y,YA,YM,Y2,XX(41),WW(52),
  4406.      *                 FV1(21),FV2(21),AA,BB,C
  4407.       EXTERNAL F
  4408.       INTEGER ISTART(4),LEN(4),KE,K1,L1,L2
  4409.       SAVE ISTART,LEN,XX,WW
  4410.       DATA ISTART/0, 7, 17, 31/,LEN/7, 10, 14, 21/
  4411.       DATA XX(  1)/.0                       /
  4412.       DATA XX(  2)/.25000000000000000000E+00/
  4413.       DATA XX(  3)/.50000000000000000000E+00/
  4414.       DATA XX(  4)/.75000000000000000000E+00/
  4415.       DATA XX(  5)/.87500000000000000000E+00/
  4416.       DATA XX(  6)/.93750000000000000000E+00/
  4417.       DATA XX(  7)/.10000000000000000000E+01/
  4418.       DATA XX(  8)/.37500000000000000000E+00/
  4419.       DATA XX(  9)/.62500000000000000000E+00/
  4420.       DATA XX( 10)/.96875000000000000000E+00/
  4421.       DATA XX( 11)/.12500000000000000000E+00/
  4422.       DATA XX( 12)/.68750000000000000000E+00/
  4423.       DATA XX( 13)/.81250000000000000000E+00/
  4424.       DATA XX( 14)/.98437500000000000000E+00/
  4425.       DATA XX( 15)/.18750000000000000000E+00/
  4426.       DATA XX( 16)/.31250000000000000000E+00/
  4427.       DATA XX( 17)/.43750000000000000000E+00/
  4428.       DATA XX( 18)/.56250000000000000000E+00/
  4429.       DATA XX( 19)/.84375000000000000000E+00/
  4430.       DATA XX( 20)/.90625000000000000000E+00/
  4431.       DATA XX( 21)/.99218750000000000000E+00/
  4432. C   NUMBER OF NODES 13
  4433.       DATA WW(1)/1.303262173284849021810473057638590518409112513421E-1/
  4434.       DATA WW(2)/2.390632866847646220320329836544615917290026806242E-1/
  4435.       DATA WW(3)/2.630626354774670227333506083741355715758124943143E-1/
  4436.       DATA WW(4)/2.186819313830574175167853094864355208948886875898E-1/
  4437.       DATA WW(5)/2.757897646642836865859601197607471574336674206700E-2/
  4438.       DATA WW(6)/1.055750100538458443365034879086669791305550493830E-1/
  4439.       DATA WW(7)/1.571194260595182254168429283636656908546309467968E-2/
  4440. C   NUMBER OF NODES 19
  4441.       DATA WW(8)/1.298751627936015783241173611320651866834051160074E-1/
  4442.       DATA WW(9)/2.249996826462523640447834514709508786970828213187E-1/
  4443.       DATA WW(15)/5.542699233295875168406783695143646338274805359780E-2/
  4444.       DATA WW(10)/1.680415725925575286319046726692683040162290325505E-1/
  4445.       DATA WW(16)/9.986735247403367525720377847755415293097913496236E-2/
  4446.       DATA WW(11)/1.415567675701225879892811622832845252125600939627E-1/
  4447.       DATA WW(12)/1.006482260551160175038684459742336605269707889822E-1/
  4448.       DATA WW(13)/2.510604860724282479058338820428989444699235030871E-2/
  4449.       DATA WW(17)/4.507523056810492466415880450799432587809828791196E-2/
  4450.       DATA WW(14)/9.402964360009747110031098328922608224934320397592E-3/
  4451. C   NUMBER OF NODES 27
  4452.       DATA WW(18)/6.300942249647773931746170540321811473310938661469E-2/
  4453.       DATA WW(28)/1.239572396231834242194189674243818619042280816640E-1/
  4454.       DATA WW(19)/1.261383225537664703012999637242003647020326905948E-1/
  4455.       DATA WW(25)/1.235837891364555000245004813294817451524633100256E-1/
  4456.       DATA WW(20)/1.273864433581028272878709981850307363453523117880E-1/
  4457.       DATA WW(26)/1.148933497158144016800199601785309838604146040215E-1/
  4458.       DATA WW(29)/2.501306413750310579525950767549691151739047969345D-2/
  4459.       DATA WW(21)/8.576500414311820514214087864326799153427368592787E-2/
  4460.       DATA WW(30)/4.915957918146130094258849161350510503556792927578E-2/
  4461.       DATA WW(22)/7.102884842310253397447305465997026228407227220665E-2/
  4462.       DATA WW(23)/5.026383572857942403759829860675892897279675661654E-2/
  4463.       DATA WW(27)/1.252575774226122633391477702593585307254527198070E-2/
  4464.       DATA WW(31)/2.259167374956474713302030584548274729936249753832E-2/
  4465.       DATA WW(24)/4.683670010609093810432609684738393586390722052124E-3/
  4466. C   NUMBER OF NODES 41
  4467.       DATA WW(32)/6.362762978782724559269342300509058175967124446839E-2/
  4468.       DATA WW(42)/1.187141856692283347609436153545356484256869129472E-1/
  4469.       DATA WW(46)/1.533126874056586959338368742803997744815413565014E-2/
  4470.       DATA WW(33)/9.950065827346794643193261975720606296171462239514E-2/
  4471.       DATA WW(47)/3.527159369750123100455704702965541866345781113903E-2/
  4472.       DATA WW(39)/8.140326425945938045967829319725797511040878579808E-2/
  4473.       DATA WW(48)/5.000556431653955124212795201196389006184693561679E-2/
  4474.       DATA WW(34)/7.048220002718565366098742295389607994441704889441E-2/
  4475.       DATA WW(49)/5.744164831179720106340717579281831675999717767532E-2/
  4476.       DATA WW(40)/6.583213447600552906273539578430361199084485578379E-2/
  4477.       DATA WW(43)/5.999947605385971985589674757013565610751028128731E-2/
  4478.       DATA WW(35)/6.512297339398335645872697307762912795346716454337E-2/
  4479.       DATA WW(44)/5.500937980198041736910257988346101839062581489820E-2/
  4480.       DATA WW(50)/1.598823797283813438301248206397233634639162043386E-2/
  4481.       DATA WW(36)/3.998229150313659724790527138690215186863915308702E-2/
  4482.       DATA WW(51)/2.635660410220884993472478832884065450876913559421E-2/
  4483.       DATA WW(37)/3.456512257080287509832054272964315588028252136044E-2/
  4484.       DATA WW(41)/2.592913726450792546064232192976262988065252032902E-2/
  4485.       DATA WW(45)/5.264422421764655969760271538981443718440340270116E-3/
  4486.       DATA WW(52)/1.196003937945541091670106760660561117114584656319E-2/
  4487.       DATA WW(38)/2.212167975884114432760321569298651047876071264944E-3/
  4488. C
  4489. C***FIRST EXECUTABLE STATEMENT QXRUL
  4490.       K=KE+1
  4491.       IS=ISTART(K)
  4492.       KS=LEN(K)
  4493.       LDL= XU-XL
  4494.       BB= LDL*0.5E0
  4495.       AA= XL+BB
  4496.       Y =0.E0
  4497.       DO 100 I=1,KS
  4498.          IF(I.GT.L1.OR.I.GT.L2)   C=BB*XX(I)
  4499.          IF(I.GT.L1)              FV1(I)=F(AA+C)
  4500.          IF(I.GT.L2)              FV2(I)=F(AA-C)
  4501. 100      Y=Y+(FV1(I)+FV2(I))*WW(IS+I)
  4502.       Y2=Y
  4503.       Y=Y*BB
  4504.       IF(L1.LT.KS)L1=KS
  4505.       IF(L2.LT.KS)L2=KS
  4506.       IF(KE.NE.K1)RETURN
  4507.       YA=0.E0
  4508.       DO 25 I=1,KS
  4509.   25    YA=YA+(ABS(FV1(I))+ABS(FV2(I)))*WW(IS+I)
  4510.       Y2=Y2*0.5E0
  4511.       YM=0.E0
  4512.       YA=YA*ABS(BB)
  4513.       DO 27 I=1,KS
  4514.    27 YM=YM+(ABS(FV1(I)-Y2)+ABS(FV2(I)-Y2))*WW(IS+I)
  4515.       YM=YM*ABS(BB)
  4516.       RETURN
  4517.       END
  4518. *DECK QXRRD
  4519.       SUBROUTINE QXRRD(F,Z,LZ,XL,XU,R,S,LR,LS)
  4520. C
  4521. C            TO REORDER THE COMPUTED FUNCTIONAL VALUES BEFORE
  4522. C            THE BISECTION OF AN INTERVAL
  4523. C           PARAMETERS
  4524. C            ON ENTRY
  4525. C              F      - REAL
  4526. C                       FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  4527. C                       FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  4528. C                       DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  4529. C
  4530. C              XL     - REAL
  4531. C                       LOWER LIMIT OF INTERVAL
  4532. C
  4533. C              XU     - REAL
  4534. C                       UPPER LIMIT OF INTERVAL
  4535. C
  4536. C              Z      - REAL
  4537. C                       VECTOR CONTAINING LZ
  4538. C                       SAVED  FUNCTIONAL VALUES
  4539. C
  4540. C              LZ     - INTEGER
  4541. C                       NUMBER OF ELEMENTS IN LZ
  4542. C
  4543. C
  4544. C            ON RETURN
  4545. C              R      - REAL
  4546. C              S        VECTORS CONTAINING LR, LS
  4547. C                       SAVED  FUNCTIONAL VALUES FOR THE NEW INTERVALS
  4548. C
  4549. C              LR     - INTEGER
  4550. C              LS       NUMBER OF ELEMENTES IN R,S RESPECTIVELY
  4551. C
  4552. C***ROUTINES CALLED  F
  4553. C
  4554.       REAL F,R,S,Z,XU,XL,DLEN,CENTR
  4555.       INTEGER LR,LS,LZ
  4556.       DIMENSION R(21),S(21),Z(21)
  4557. C***FIRST EXECUTABLE STATEMENT QXRRD
  4558.       DLEN= (XU-XL)*0.5E0
  4559.       CENTR= XL+DLEN
  4560.       R(1)=  Z(3)
  4561.       R(2)=  Z(9)
  4562.       R(3)=  Z(4)
  4563.       R(4)=  Z(5)
  4564.       R(5)=  Z(6)
  4565.       R(6)=  Z(10)
  4566.       R(7)=  Z(7)
  4567.       S(1)=  Z(3)
  4568.       S(2)=  Z(8)
  4569.       S(3)=  Z(2)
  4570.       S(7)=  Z(1)
  4571.       IF(LZ.GT.11)GOTO 2
  4572.       R(8)=  F(CENTR+DLEN*0.37500000E0)
  4573.       R(9)=  F(CENTR+DLEN*0.62500000E0)
  4574.       R(10)=  F(CENTR+DLEN*0.96875000E0)
  4575.       LR=  10
  4576.       IF(LZ.NE.11)S(4)=  F(CENTR-DLEN*0.75000000E0)
  4577.       IF(LZ.EQ.11)S(4)=  Z(11)
  4578.       S(5)=  F(CENTR-DLEN*0.87500000E0)
  4579.       S(6)=  F(CENTR-DLEN*0.93750000E0)
  4580.       S(8)=  F(CENTR-DLEN*0.37500000E0)
  4581.       S(9)=  F(CENTR-DLEN*0.62500000E0)
  4582.       S(10)=  F(CENTR-DLEN*0.96875000E0)
  4583.       LS=  10
  4584.       GOTO 3000
  4585. 2     R(8)= Z(12)
  4586.       R(9)= Z(13)
  4587.       R(10)= Z(14)
  4588.       LR=  10
  4589.       S(4)= Z(11)
  4590.       S(5)= F(CENTR-DLEN*0.87500000E0)
  4591.       S(6)= F(CENTR-DLEN*0.93750000E0)
  4592.       IF(LZ.GT.14)GOTO3
  4593.       S(8)= F(CENTR-DLEN*0.37500000E0)
  4594.       S(9)= F(CENTR-DLEN*0.62500000E0)
  4595.       S(10)= F(CENTR-DLEN*0.96875000E0)
  4596.       LS=  10
  4597.       GOTO 3000
  4598. 3     R(11)= Z(18)
  4599.       R(12)= Z(19)
  4600.       R(13)= Z(20)
  4601.       R(14)= Z(21)
  4602.       LR=  14
  4603.       S(8)= Z(16)
  4604.       S(9)= Z(15)
  4605.       S(10)= F(CENTR-DLEN*0.96875000E0)
  4606.       S(11)= Z(17)
  4607.       LS=  11
  4608. 3000  RETURN
  4609.       END
  4610. *DECK QXCPY
  4611.       SUBROUTINE QXCPY(A,B,L)
  4612. C
  4613. C  TO COPY THE REAL VECTOR B OF LENGTH L   I N T O
  4614. C          THE REAL VECTOR A OF LENGTH L
  4615. C
  4616. C***REMARK  THIS ROUTINE CAN BE IMPROVED, BY CODING IT IN THE
  4617. C           ASSEMBLER LANGUAGE OF THE USED MACHINE
  4618. C***ROUTINES CALLED  (NONE)
  4619. C
  4620.       INTEGER L
  4621.       REAL A(L),B(L)
  4622. C***FIRST EXECUTABLE STATEMENT QXCPY
  4623.       DO 1 I=1,L
  4624. 1     A(I)=B(I)
  4625.       RETURN
  4626.       END
  4627. *DECK MAIN *** TEST PROGRAM ***
  4628. C  THIS DRIVER CALLS VARIOUS INTEGRATORS WITH  VARIOUS RELATIVE
  4629. C  TOLERACES TO INTEGRATE THE FUNCTION
  4630. C     F(X) = SQRT(X/(1-X)) * LOG(X)   OVER (0, 1)
  4631. C
  4632. C THE EXACT INTEGRAL IS                           -0.606789763508705...
  4633. C
  4634.       REAL F,A,B,ATOL,RTOL,TEE,Y,WORK(4600)
  4635.       INTEGER NCALLS,IER,NEVAL,LAST,IWORK(300),KKK,KEY,LIMIT,LENW,LENIW
  4636.       EXTERNAL F
  4637.       COMMON /CTEST/NCALLS
  4638.       DATA LIMIT/100/,LENIW/300/,LENW/4600/
  4639.       WRITE(6,600)
  4640. 600   FORMAT(32H1 COMPUTATION OF THE INTEGRAL OF,
  4641.      *      /38H       F(X) = SQRT(X/(1-X)) * LOG(X),
  4642.      *      /13H  OVER (0, 1),
  4643.      *      /23H  THE EXACT INTEGRAL IS,30X,
  4644.      *       22H -0.606789763508705...)
  4645.       A=0.0
  4646.       B=1.0
  4647.       WRITE(6,61)
  4648.  61   FORMAT(//22H TEST OF QXGS                //
  4649.      *  53H RELATIVE TOLERANCE FUNCTION CALLS ERROR ESTIMATE IER,
  4650.      *  13H       RESULT )
  4651.       DO 1 I=1,4
  4652.          ATOL=1E-8
  4653.          RTOL=10E0**(-I)
  4654.          NCALLS=0
  4655.          CALL  QXGS (F,A,B,ATOL,RTOL,Y,TEE,IER,LIMIT,
  4656.      *               LENIW,LENW,LAST,IWORK,WORK)
  4657.          WRITE(6,60)I,NCALLS,TEE,IER,Y
  4658. 60       FORMAT(4X,6H10**(-,I2,1H),I17,9X,E7.2,4X,I2,F20.8)
  4659.  1    CONTINUE
  4660.       DO 2 KKK = 1,5
  4661.        KEY = KKK - 1
  4662.        WRITE(6,62) KEY
  4663.  62    FORMAT(//23H TEST OF  QXG , KEY =    ,I4   //
  4664.      *  53H RELATIVE TOLERANCE FUNCTION CALLS ERROR ESTIMATE IER,
  4665.      *  13H       RESULT )
  4666.       DO 2 I=1,4
  4667.          ATOL=1E-8
  4668.          RTOL=10E0**(-I)
  4669.          NCALLS=0
  4670.          CALL  QXG  (F,A,B,ATOL,RTOL,KEY,Y,TEE,IER,LIMIT,
  4671.      *               LENIW,LENW,LAST,IWORK,WORK)
  4672.          WRITE(6,60)I,NCALLS,TEE,IER,Y
  4673.  2    CONTINUE
  4674.       STOP
  4675.       END
  4676. *DECK F
  4677.       REAL FUNCTION F(X)
  4678.       REAL X
  4679. C        INTEGRAND FUNCTION, IN THE INTERVAL (0,1) THE EXACT INTEGRAL
  4680. C   TO 15 DIGITS IS
  4681. C          -0.60678 97635 08705D0
  4682. C
  4683.       COMMON /CTEST/NCALLS
  4684.       NCALLS=NCALLS+1
  4685.       IF (X.EQ.0.0) GO TO  10
  4686.       IF (X.EQ.1.0) GO TO  10
  4687.       F=SQRT(X/(1-X))*LOG(X)
  4688.       RETURN
  4689. 10    F=0.0
  4690.       RETURN
  4691.       END
  4692.